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ABSTRACT 


Modifications  and  improvements  to  the  Nabetani-Rankin 
method  of  sequential  layer  addition  for  a  one-dimensional 
earth  are  presented  and  applied  to  several  representative 
models  to  demonstrate  the  technique  developed  for  inversion 
of  magnetotelluric  data.  The  inversion  process  was  aided 
by  the  use  of  computer  graphics  and  interactive  programming. 
It  was  found  that  there  are  general  rules  to  be  followed 
for  choosing  layer  parameters  and  once  these  methods  are 
learned,  the  advantage  and  speed  of  interactive  program¬ 
ming  can  be  utilized.  This  method  is  inherently  stable 
compared  with  the  "least-squares"  method  of  curve  fitting. 
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CHAPTER  1 

INTRODUCTION 

1.1.  Purpose  of  the  Investigation 

The  generalized  magnetotelluric  (MT)  method  consists  of 
measuring  the  horizontal  components  of  the  naturally-occur¬ 
ring  electromagnetic  field  variations  at  the  surface  of  the 
earth  and  interpreting  the  sub-surface  structure  in  terms  of 
the  apparent  resistivity  or  the  surface  impedance.  Direct 
interpretation  of  MT  data,  consisting  of  master-curve  match¬ 
ing  by  trial  and  error,  can  be  very  tedious,  especially  if 
the  number  of  parameters  involved  is  large.  Due  to  the 
availability  of  fast  electronic  computers  and  reliable  record 
ing  equipment,  the  possibility  exists  of  automatic  inversion 
of  MT  data,  which  in  turn  may  be  a  more  reliable  method  and 
can  save  considerable  computation  time.  Existing  data 
inversion  techniques  have  been  based  on  "least-squares"  meth¬ 
ods  of  curve  fitting,  in  which  the  mathematical  and  computa¬ 
tional  procedures  tend  to  obscure  the  physical  basis  of  the 
method  and  can  lead  to  ambiguous  or  faulty  results.  The 
method  of  data  inversion  presented  in  this  work  makes  use  of 
the  exact  magnetotelluric  expressions,  which  are  intimately 
connected  with  the  physical  model  and  can  thus  match  as 
closely  as  the  measured  data  warrant. 

1.2.  Historical  Review  of  the  Magnetotelluric  Method 

Tikhonov  (1950),  Rikitake  (1950,  1951)  and  Kato  and 
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Kikuchi  (1950)  observed  the  existence  of  a  correlation 
between  the  geomagnetic  and  telluric  field  variations  at  the 
surface  of  the  earth,  and  suggested  that  it  would  be  useful 
for  depth  sounding.  The  method  did  not  become  popular  until 
Cagniard's  definitive  paper  (1953),  which  presented  a  graph¬ 
ical  technique  for  interpreting  magnetotelluric  field  data 
under  the  assumption  that  the  earth  is  horizontally  strati¬ 
fied  (that  is,  one-dimensional)  and  that  the  source  fields 
are  plane  em  waves.  Wait  (1954)  and  Price  (1962)  considered 
the  problem  of  finite  source  dimensions  and  showed  that 
Cagniard's  results  are  valid  only  if  the  fields  do  not  vary 
significantly  over  a  horizontal  distance  of  the  order  of  a 
skin  depth.  However,  Madden  and  Nelson  (1964),  considering 
realistic  earth  models,  concluded  that  Cagniard's  plane  wave 
assumption  is  valid  in  most  cases. 

Magnetotelluric  measurements  have  been  made  by  many 
investigators  (Cantwell,  1960;  Srivistava  et  al.,  1963; 

Vozoff  et  al.,  1963;  Tikhonov  and  Berdichevskii,  1966;  Swift, 
1967;  Peeples,  1969;  Reddy  and  Rankin,  1971;  and  others). 

When  the  earth  is  anisotropic  or  inhomogeneous,  the 
resistivity  must  be  treated  as  a  tensor  quantity.  Chetaev 
(1960)  ,  Cantwell  (1960)  ,  Kovtun  (1961)  ,  Bostick  and  Smith 
(1962),  Swift  (1967),  Morrison  et  al.  (1968)  and  others  have 
shown  how  to  compute  the  tensor  components  using  spectral  and 
statistical  techniques. 

The  nature  of  the  MT  field  for  cylindrical  inhomogeneit¬ 
ies  has  been  studied  by  several  investigators,  representative 
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works  being  those  of  Neves  (1957),  who  studied  the  response  of 
the  MT  field  to  dipping  interfaces  using  finite  difference 
techniques;  d'Erceville  and  Kunetz  (1962),  who  solved  analyt¬ 
ically  for  the  vertical  contact  fault;  and  Rankin  (1962) ,  who 
solved  analytically  for  the  vertical  dike. 

Swift  (1967),  following  the  approach  of  Madden  (1966), 
solved  for  two-dimensional  problems  by  use  of  a  transmission 
line  analogy  and  used  the  method  to  study  an  electrical 
conductivity  anomaly  in  the  southwest  United  States. 

There  have  been  many  recent  attempts  at  inverting  MT 
data  (for  example,  Chetaev,  1966).  The  most  common  method  is 
that  of  "least-squares"  fitting  of  experimental  curves  (see 
Marquardt,  1963) .  Least-squares  approaches  have  been  adopted 
by  Nelder  and  Mead  (1964),  Tikhonov  (1965),  Wu  (1968), 

Patrick  (1969) ,  Patrick  and  Bostick  (1969) ,  Dmitriev  (1970) , 
Laird  and  Bostick  (1970) ,  and  Word  (1970)  . 

In  this  work,  the  inversion  technique  of  Nabetani  and 
Rankin  (1969)  is  used;  the  theoretical  development  and  applic¬ 
ation  follow  in  subsequent  chapters. 

1.3.  Existing  Methods  of  Inverse  Analysis 

Tikhonov  (1965)  showed  the  uniqueness  of  the  solution  of 
the  inverse  problem  of  MT  sounding.  That  is,  if  the  impedance 
which  in  general  can  be  written 

Z  =  Z  ( p ;  oj)  (1.1) 


or  specifically 
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Z  =  z(oj) 

(1.2) 

then 

P  =  P  (z) 

(1.3) 

is  a  unique  dependance  of  resistivity  p  on  depth  z. 
Thus , 


Z  (a))  =  A  [ p  (  z  )  ] 


(1.4) 


where  A  is  a  nonlinear  operator. 

Inverse  methods  that  have  employed  the  "least-squares" 
technique  involve  the  computation  of  either  apparent  resisti¬ 
vity  or  surface  impedance.  In  this  section,  the  theory  of 
least-squares  fitting  will  be  developed  briefly. 

Letting  p  .  be  the  calculated  and  p  .  be  the  observed 
ci  oi 

apparent  resistivity  at  period  T^,  where  the  apparent  resisti¬ 
vity  is  a  function  of  p,z,T^  (which  will  be  developed  explic¬ 
itly  in  Chapter  2)  , one  computes  the  parameters  that  minimize 
the  function 


z 

i 


(1.5) 


where 


V  -  ^(p.,h.) 


(1.6) 


with  p being  the  resistivity  and  a  thickness  parameter 
associated  with  the  j ' th  value  of  p.  Local  minima  occur 
where 


3Y/3Xj  =  0  ,  j  =  1,2,  ...  ,  2n-l  (1.7) 


J 
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where  is  the  j ' th  parameter  of  y.  That  is, 


A1  pl'  '  An  pn'  An+1  "  hl'  '  A 


,  =  h  .  (1.8) 

zn-l  n- I 


Writing 


pci  f(pl'  '  pn '  hl'  •••  '  hn-l '  Ti)  "  f(A;T)  (1’9) 


if  f  is  linear  in  the  A's  and  the  A's  are  linearly-independent , 
the  ¥  surfaces  are  ellipsoids: 


¥  - 


(pl  -  ap  (p2  -  a2)2 

“ ““ ““ — “ “ ~ “  "I"  -  "(■ 


.  + 


(hn-l  -  a2n-l> 


2n-l 


(1.10) 


If  f  is  not  linear,  the  surfaces  are  distorted;  however,  in 
the  immediate  vicinity  of  the  minimum  ¥ ,  the  function  f(X;T) 
behaves  in  a  linear  manner. 

For  a  change  in  one  of  the  parameters,  there  will  be  a 
corresponding  change  in  p  Consider  the  finite  difference 

approximation 


9pci/9pj  ~  Apci/APj  =  1/Apj  f(fV  •••  '  Pj-1'  Pj  +  Apj' 


Pj  +  1'  • • •  *  pn'  hi'  • * •  '  hn-l;  Ti} 


(1.11) 


and  similarly  for  changes  in  the  h^ . 
To  first  order, 


pci> 


9  p  .  /  9  p  dp,  +  ...  + 
ci  1  1 


9p  ./9h  dh 
ci  n  n 


(1.12) 


. 


' 


■ 


6 


and  thus 


*  =  £  (dp.)2 

i  cl 


(1.13) 


Defining  the  matrices 


D  =  (dp1, 


'  dhn> 


(1.14) 


P  = 


9p 


cl 


9  p 


9p 


cn 


9  p . 


3p 


cl 


9h 


3p 


cn 


9h 


n 


(1.15) 


R  (pol  P^1 '  *  *  *  '  p 


cl 


on 


P  _  _ ) 
cn 


(1.16) 


One  can  then  write  (1.12)  as 


T  T 

P  P  D  =  P1  R 


(1.17) 


And  solving  for  D  by  inversion, 


T  —  T 

D  =  (P  P)  !P  R 


(1.18) 


However,  since  p  .  is  nonlinear  in  p  and  h,  neglecting  high- 

C 1 

er  order  derivatives  in  (1.12)  may  not  lead  to  convergence  if 
the  norm  of  D  is  large.  This  can  be  seen  by  considering 
(1.17),  in  which  a  large  value  of  D  would  require  a  small 
value  of  the  norm  of  P.  P  is  composed  of  terms  resembling 
(1.11),  which  would  itself  be  quite  sensitive  to  changes  in 


the  values  of  its  parameters,  since  its  coefficients  must  of 
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necessity  be  very  small.  Thus,  the  operation  (1.18)  is  unre¬ 
liable  and  the  process  beset  by  inherent  instability. 

Most  algorithms  for  least-squares  estimation  of  nonlinear 
parameters  have  been  based  on  either  a  Taylor  series  expan¬ 
sion  of  (1.9)  or  modifications  of  the  method  of  steepest  desc¬ 
ent. 

1.3. a  Expansion  of  f(X;T)  in  a  Taylor  Series 

In  this  method,  corrections  to  the  parameters  are  calcu¬ 
lated  on  the  assumption  of  local  linearity  of  the  function 
about  the  point  A  =  AQ  +  6A: 


(1.19) 


and  truncating  after  the  last  linear  term  (since  we  assume 
that  f  is  linear  in  X  and  thus  the  higher  derivatives  are 
zero),  (1.19)  can  be  written 


(1.20) 


f 


where  f  is  the  original  (uncorrected)  function.  Now  5A  ap¬ 


pears  linearly  in  (1.19)  and  (1.20)  and  can  be  found  by  set¬ 
ting 


3 y/3  (6 A ^ ) 


0 


(1.21) 


for  all  j. 


That  is,  by  solving 


A  6A 


g 


(1.22) 


'  * 
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where 


A 


P 


g 


(1.23) 

(1.24) 


(1.25) 


The  corrections  6  A  in  each  iteration  must  be  small;  otherwise 
the  extrapolation  may  be  beyond  the  region  where  f  can  be 
considered  linear  and  the  method  may  not  lead  to  convergence. 
The  conditions  (1.21)  lead  to  a  stationary  value  of  Y ,  which 
may  not  be  a  minimum  for  all  parameters  simultaneously.  One 
most  important  requirement  is  a  good  "first  guess"  for  the  A_. 
and  even  so,  convergence  is  slow. 

1.3.b.  Another  technique  of  solution  is  the  modified  Newton- 
Raphson  method,  where  the  second-order  term  in  (1.19)  is 
retained : 


(1.26) 


It  has  been  found,  though,  that  this  method  does  not  always 
decrease  ¥  at  each  iteration.  (See  Laird  and  Bostick,  P.  21.) 

1.3.c.  Gradient  Methods 

One  moves  in  the  direction  of  the  negative  gradient  of  Y 
(that  is,  in  the  direction  of  minimum  ¥) .  The  step  taken  in 
each  iteration  is  thus 


6  A 


[ 3^/3  A  x 


(1.27) 


1 

I H  1 1 1  I  ■ 
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The  control  of  step  direction  and  size  is  important  here,  for 
a  bad  choice  could  lead  to  divergence.  Convergence  in  the 
neighbourhood  of  minimum  ¥  is  as  a  rule  extremely  slow. 
Marquardt  (1963)  has  devised  a  method  for  determination  of 
step  size  and  direction  simultaneously,  but  uniqueness  is  not 
guaranteed.  (See  Wu,  1968.) 

1.4.  Outline  of  Thesis 

In  Chapter  2,  the  basic  equations  of  the  magnetotelluric 
field  are  derived  and  applied  to  a  layered  half-space.  The 
Nabetani-Rankin  inverse  method  is  introduced  and  the  relation 
ships  to  be  used  in  applying  it  are  derived. 

The  method  of  sequential-layering  and  the  theory  behind 
the  practical  application  are  presented  in  Chapter  3. 

In  Chapter  4,  computer  graphics  diagrams  of  curve  match¬ 
ing  by  the  method  of  sequential-layering  are  shown.  Several 
representative  earth  models  are  presented  to  offer  a  feeling 
for  the  method.  Chapter  4  is  followed  by  conclusions  and 
suggestions  for  further  work,  references  and  three  Appendices 


* 


10 


CHAPTER  2 

MAGNETOTELLURIC  THEORY 

In  the  em  system  of  units,  Maxwell's  equations  can  be 
written 


V 

X 

E  = 

-> 

-9H/9t 

(2.1) 

V 

X 

H  = 

-> 

4ttj  +  9  D/9 1 

(2.2) 

V 

• 

E  = 

47Td 

(2.3) 

V 

• 

H  = 

0 

(2.4) 

where  J  is  the  electric  current  density  and  is  related  to  the 
electric  field  by 

j  =  qe  =  E/p  (2.5) 

where  o  is  the  electric  conductivity,  p  the  electric  resisti¬ 
vity,  and  d  the  charge  density  of  the  medium.  One  can  consi¬ 
der  a  single  Fourier  component  of  a  general  em  disturbance, 
writing  the  time-dependent  fields  as 

{  E  }  =  {  E°  }eiWt  (2.6) 

B  H0 

In  equation  (2.2),  applied  to  media  such  as  the  earth,  dis¬ 
placement  currents  are  negligible  in  comparison  with  conduc¬ 
tion  currents  in  the  frequency  range  of  interest  (cq<1045£^_l)  . 

Sec . 

Since  we  are  considering  media  of  finite  conductivity,  local 
concentrations  of  charge  can  also  be  neglected.  Hence,  the 
previous  equations  may  be  rewritten  as: 


' 


. 
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v  x  £  =  -aS/at  (2.7) 

V  x  H  =  4ttJ  (2.8) 

V  •  £  =  0  (2.9) 

V  *  H  =  0  (2.10) 


To  develop  the  MT  equations,  one  solves  Maxwell's  equa¬ 
tions  subject  to  the  following  assumptions.  Plane  em  waves 
impinge  perpendicular  to  the  surface  of  the  earth  and  propa¬ 
gate  vertically  downward.  This  latter  point  is  true  for  a 
non-zero  angle  of  incidence,  due  to  the  relatively  large 
refractive  index  of  the  earth  with  respect  to  the  air  for 
electromagnetic  radiation.  The  co-ordinate  system  used  is 
right-handed  Cartesian  with  z  vertically  downward  and  the  x-y 
plane  on  the  surface  of  the  earth.  Since  we  are  not  concerned 
with  magnetic  materials,  we  can  consider  the  magnetic  permea¬ 
bility  to  be  approximately  that  of  free  space.  That  is, 


Assuming  the  em  fields  polarized  such  that 


—r 

E  = 

(E  ,0,0) 

A 

(2.11) 

-¥ 

H  = 

(0,H  ,0) 

(2.12) 

and  making  use  of  the 

vector 

identity 

V  X  (V 

X  E)  = 

V  (V  •  E)  -  V2  E 

(2.13) 

From  (2.7)  and  (2.6), 

V  X  E  =  -3H/3t 


-iwH 


(2.14) 


ai  is rfT 


' 

' 
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From  (2.8)  and  (2.5), 

V  X  H  =  4ttJ  =  4ttqE  (2.15) 

V  X  (V  X  2)  =  -iu)  (V  X  S) 

=  -4iraia)£  (2.16) 

Making  use  of  (2.9),  this  equation  can  be  written 

V2  2  =  47Tai^E  (2.17) 

with  an  identical  equation  for  H.  The  propagation  equations 

for  electromagnetic  waves  in  the  earth  are  thus 

-> 

E 

(V2  -  k2)  (+}  =  0  (2.18) 

H 

where  k2  =  4Traiu)  (2.19) 

Solving  equation  (2.18),  for  example,  for  the  electric  field 
for  the  special  case  of  an  isotropic  half-space,  where  the 
partial  derivatives  with  respect  to  x  and  y  are  zero  because 
of  symmetry, 

V2  E  =  d2E  /dz2  (2.20) 

A 

which  has  the  general  solution 

E  =  Aekz  +  Be_kz  (2.21) 

X 

The  first  term  corresponds  to  an  upward-travelling  and  the 
second  term  to  a  downward-travelling  wave.  In  this  partic¬ 
ular  case,  there  is  no  reflected  (upward-travelling)  wave; 


. 


-'oiJifiq  airfi  nl  .3vsw  pnxXlsveia-b^Bwnwob  b  oj  im»u  mm  A 
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indeed,  for  the  em  fields  to  remain  finite  for  large  values 
of  z,  the  waves  must  decay  with  depth.  Hence,  A  =  0.  There¬ 
fore  , 


E 

x 


(2.22) 


and  from  (2.14) 


io)  dz 


k  _  -kz 
- —  Be 
10) 


(2.23) 


One  solves  for  p  by  using  equation  (2.19)  and  (2.23), 


=  2T 


x 


H 


(2.24) 


and  the  depth  at  which  an  electromagnetic  Fourier  component 
has  diminished  in  magnitude  to  1/e  of  its  original  value 
(known  as  the  "skin  depth") 

&  =  27  (PT) 1/2  (2.25) 

In  the  practical  system  of  units,  one  measures  the  mag¬ 
netic  field  in  gammas,  the  electric  field  in  millivolts/km. , 
the  period  of  wave  oscillation  in  seconds,  and  the  resistiv¬ 
ity  in  ohm-metres. 
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1  gamma 


_  5 

10  em  cgs 


1  mv/km 


1  em  cgs 

5 

10  em  cgs 

i  n1  1 

10  em  cgs 


1  km 


1  ohm-M 


In  this  system  of  units 


P 


0.2  T 


(2.26) 


H 

y 


6 


(2.27) 


The  case  of  a  finite  source  (non-plane  waves)  has  been 
dealt  with  by  Wait  (1954)  and  more  completely  by  Price  (1962) 
Price's  analysis  can  be  found  in  APPENDIX  A.  Madden  and 
Nelson  (1964)  have  shown,  however,  that  the  plane  wave  assump 
tion  is  valid  for  the  frequencies  of  interest  in  the  MT 
method. 

Cagniard  (1953)  developed  relationships  for  the  apparent 
resistivity  in  analogy  to  equation  (2.26)  for  the  case  of 
two-  and  three-layer  earth  models  and  outlined  a  method  for 
interpreting  experimental  results  by  matching  with  master 
curves.  Master  curves  for  2,3,4,  and  5-layer  earth  models 
have  been  computed  by  various  authors  (such  as  Yungul,  1961; 
Hasegawa,  1962;  and  Srivistava,  1967).  However,  the  labori¬ 
ous  algebraic  work  involved  can  be  avoided  by  a  simple  itera¬ 
tive  process  for  a  general  n-layered  model  using  the  digital 
computer  directly.  Adopting  the  co-ordinate  system  chosen 
earlier,  so  that  the  em  fields  are  given  by  (2.11)  and  (2.12) 
and  the  notation  for  the  earth  cross-section  in  FIGURE  1,  the 
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FIGURE  1: 


n  +  1  layer  earth  model 
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wave  equation  for  the  electric  field  in  the  m'th  layer  can  be 
written  (from  2.18) 


d2E 


x,m  _  ,2 


dz 


=  k  E 

m  x,m 


where 


k2  =  4 7T a  iw  =  4TTia)/p 

m  m  m 


which  has  the  general  solution 


x,m 


A  e 
m 


-k  (z  -  z  ->  ) 
m  m-1 


+  B  e 
m 


k  ( z  -  z  i  ) 
m  m-1 


The  magnetic  field,  given  by  (2.14),  is 


H 


y,m 


,  dE 

_ 1  /  x,mx 

i  w  ldz  ; 


k  -k  (z  -  z  ,  ) 

J{Ae  m  m-X 

iw  m 


k  (z  -  z  ,  ) 

B  e  m  m-1  ] 

m  J 


Applying  the  boundary  conditions,  that  the  tangential 
ents  of  E  and  H  are  continuous  across  a  discontinuity, 


x,m 


::=z 


m 


=  E 


x,m+l  _ 


z=z 


m 


H 


y,m 


z=z 


m 


=  H 


y  ,m+l 


z=z 


m 


hence 


-k  h  k  h 

Ae  mm  +  Bemm  =  A  +  B 

m  m  m+ 1  m+ 1 


(2.28) 

(2.29) 

(2.30) 


(2.31) 

compon- 

(2.32) 

(2.33) 


(2.34) 
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-k  h  k  h 

Ae  mm  -  Bemm 
m  m 


m+1  r  A  , 

k  iAm+l  Bm+1} 

m 


(2.35) 


h  =  z  -  z 
m  m  m-1 


(2.36) 


hi  =  zi 


Solving  for  the  coefficients, 


k  h 
e  m  m 

-  {  (k  +  k  ..  )  A  +  (k  +  k  ,.)B  (2.37) 

m  2k  m  m+1  m+1  m  m+1  m+1 

m 


-k  h 
m  m 


B 


m  2k 


m 


{  (k  +  k  ,  n  )  A  , 
1  m  m+1  m+1 


(k  +  k  )  B  } 
m  m+1  m+1 ; 


(2.38) 


For  solutions  to  remain  finite  as  ,  we  demand  that  B  .  =0 

n+1 

in  (2.30)  (for  the  same  reason  as  given  following  equation 
(2.21)).  If  one  puts  An+^  =  1'  since  in  any  event  it  cancels 
out  in  the  final  ratios, 


A 

n 


k  + 
n+1 


2k 


k  k  h 
n  n  n 
-  e 


n 


B 

n 


k  - 
n 


n+1 


-k  h 
n  n 


2k 


n 


At  the  surface  of  the  earth, 


E 

x 


z=0 


+  B 


1 


(2.39) 


(2.40) 


(2.41) 

(2.42) 


In  analogy  with  equation  (2.26),  one  defines  the  apparent 


; 


.  . 
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resistivity  p^: 


pa  =  (0.2  T) 


x 


H 


(2.43) 


z=0 


and  thus 


p  =  (0.2  T) 


Ai +  Bi 


is  (Ai  -  Bi> 


(2.44) 


THE  INVERSE  METHOD 


For  our  purposes,  it  is  more  convenient  to  deal  with  the 


intrinsic  impedance 


Z(a>)  = 


x 


lwH 


(2.45) 


z=0 


which  for  an  isotropic  half-space  from  equation  (2.23)  can  be 


written 


Z(u>) 


1 

k 


1 


(47Tai(jo) 


1/2 


(2.46) 


or 


Z(T) 


JL 

2  77 


(pT/2) 


1/2  -iir/4 
'  e 


(2.47) 


It  should  be  noted  in  this  expression  that  the  electric  field 
"leads"  the  magnetic  field  by  a  phase  angle  of  tt/4  radians. 
For  a  layered  earth,  one  can  write  from  (2.41)  and  (2.42), 


. 
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FIGURE  2: 


The  matrix  equation  (2.52) 
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Z  (oo)  = 


A0  +  Bo 


o  ^Ao  V 


Ai  +  Bi 

ki  (Ai  "  Bi} 


(2.48) 


since  from  (2.34)  and  (2.35)  at  the  earth  surface  (z=0) , 


Ao  +  Bo  -  Ai  +  B1 


(2.49) 


ko  (Ao  -  V  =  kl  (A1  '  V 


(2.50) 


where  is  the  propagation  constant  for  free  space  (k^^/c)  . 
The  ( 2n+2 )  equations  (2.37),  (2.38),  (2.39),  (2.40),  (2.49), 

and  (2.50)  can  be  written  in  matrix  notation 


M-(x)=y  (2.51) 

or  (see  APPENDIX  B)  as  in  (2.52),  FIGURE  2.  The  solution  of 
(2.51)  is,  by  inversion, 

x  =  M-1  (y)  (2.53) 


By  the  use  of  determinants,  one  writes 


x 


1 


det  M  i  i 
det  M  y 


(2.54) 


where  Mi i  is  the  appropriate  cofactor.  Or,  from  (2.52), 


-  (B0/Aq) 


det  Mi i 

det  M 


(2.55) 


Therefore , 


x 


100H 


z=0 


1  det  M  -  det  Mi i 
kg  det  M  +  det  Mi i 


(2.56) 


21 


Because  of  the  diagonal  nature  of  the  matrix  M,  the  solution 
for  the  required  determinants  is  simplified  (See  APPENDIX  B) . 
For  example,  for  a  2-layered  earth  (n=l) , 


x 


lWH 


z=0 


11+  Kie 


-2kizi 


-2k1z 

1  -  k  !  e  1  1 


where 


k.  -  k.._ 
3  3+1 

k  .  +  k .  _ 
1  3  +  1 


For  a  3-layered  earth  (n=2) , 


(2.57) 


(2.58) 


x 


iwH 


z=0 


=_1  l+K,e'2klZl+K,e~2klZl~2k2(z2~Zl)+K1K?e'2k2(22'Zl) 

kl  l-K1e-2k‘z‘-(c2e-2k*z>-2k2(z*-z>)+KlKle-2k2(Zi-z‘) 

(2.59) 


One  can  define 


det  M  -  det  Mi i 
n  det  M  +  det  Mj i 


A  o  +  Bp 

A  o  —  Bo 


(2.60) 


where 


K°  +  K1  +  ...  + 

_ _n _ n _ n 

K°  -  K1  +  -  ...  +  (-1) nK^ 
n  n  n 


n  i 
1=0  n 

n  i  i 

Z  [(-1)  IT] 


(2.61) 


for  an  n-layered  earth;  and 
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K 0  =  1 

n 


K1 

n 


n  p 

Z  k  exp  [-2  Z  k  h  ] 
1  =  1  1  1 


P=i  p 


n-1  n  q  p 

2  £  kk  exp  [-2  Zkh  +  2  Z  k  h  ] 

r.  rr  ,  ,  m  III 


k!  =  p=.l  q=l  P  q 


n 


1  =  1  1  1 


m=l 


Km  = 
n 


n+l-m  n+2-m 
Z  Z.  . 

Px=l  P0=P-,+l 


n 

Z 


2  *m-l 


Pm=Pm-i+1  3=1  Pm  t  =  1 


m  m-X  m-l+x  pi 

II  (k„  )n  (exp  [2  (-1)  Z  k  h  ]) 


qi=Xqi  qx 


m 

X  exp [-2  Zkh] 
q  =1  q  q 


n  n  n/2 

K  =  n  (k  )  exp [-2  Z  (k9  h_  )] 
j=l  J  1*1 


,  n  even 


n 


(n-1)  / 2 


=  II  (k.)  exp  [-2  Z  (k  _  n  h  ]  ,  n  odd 

j=l  J  i=l  1  Zl  1 


If  one  defines  (See  2.58) 


(2.62) 


K  . 

D 


✓oT 


/57 

J 


+  /o 


j+1 


/p 


i±i. 


/p 


i/vT  +  ■'pj 


(2.63) 


and 


H  =  47Th1//p1 


(2.64) 


j  /  Px/Pj 


W.  = 
D 


hj/hl 


(2.65) 
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Then, 


where 


-2kjhj  =  -/2i/T  H  Wj 

/2i  =  /2  e^/4  -  1  +  i 


(2.66) 


Thus,  if 


Q  =  /2i/T  H 


(2.67) 


then 


-2k . h .  =  -Q  w . 

3D  3 


(2.68) 


Equation  (2.62)  then  becomes 


K 


m 

n 


n-m+1  n-m+2 
=  I  E 

Pi=i  p2=px+i 


n 

Z 

p  =p  , +1 

m  *m-l 


m 

n  (k  ) 

.  ,  p 

3=1 


m-1  .  p  p  Q 

X  {  n  (exp [2 (-lT  5  1  (W  ) ] ) exp [-2  Im  (v»  )]} 

?=1  q?=l  qC  qm=l  qm 

( 2.69) 


If  we  further  define  (and  use  2.60) 


1  -  r 


V  = 


n 


1  - 


An  +  B 


An  -  B 


i  +  r 

n 


s_  =  _ 


1  + 


A  o  +  Bp 

Ao  —  Bq 


(B0/Aq ) 


(2.70) 


From  (2.55) 


det  ! 


V 


det  M 


(2.71) 


. 


From  (2.70)  and  (2.61), 
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n 


1  - 


V  = 


1  + 


Z  K- 
i=0 
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n 

Z 

(-1)1 

K1 
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o 

II 

•H 

n 

E 

K1 
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ll 

-H 

n 

n 

E 

(-1)1 

K1 

o 

ll 

-H 

n 

n  .  . 

Z  [  (-1)1  -  1  ]  K* 

i=0 _ _ 

z  [  (-1)1  +  i  ]  K^; 

i=0 


For  n=l,  a  2-layer  earth  (from  2.61  and  2.60), 


(2.72) 


K? 


+  K: 


kT 


k} 


-  kj 


i  + 

1  3  kT 


(2.73) 


1  +  Bq/Aq 

1  -  Bo/Ao 


(2.74) 


Thus ,  from  (2.70) , 

V  =  -B„/A0  =  -k! 


=  -k  i  e 


-QWi 


(2.75) 


For  n=2 , 


1  +  K\  +  K|_ 

I  1  ici  4‘ 


+  K2 


1  + 


1  - 


K- 


+ 

’Ki" 


s 


1  -  Ki 


(2.76) 


so  that  from  (2.7ft) , 


V  =  - 


1  + 


kJ 


-QW  i  , 

Kie  +  K2e 


~Q (W i  +  W2) 


1  +  k  i  <2e 


-QW2 


(2.77) 

(2.73-) 
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and  the  "V-function"  for  any  value  of  n  can  be  obtained  from 
(2.73)  and  (2.62).  It  can  be  regarded  as  the  negative  of  the 
ratio  of  total  reflected  to  total  transmitted  radiation  at 
the  surface  of  the  earth  (z=0) .  It  should  be  noted  that  for 
a 2<o i,  k2<ki  and  both  and  are  greater  than  zero. 

Thus,  V  =  — B o /A o  and  the  function  will  be  negative.  That  is, 
it  will  behave  like  a  curve  of  -(apparent  resistivity).  Thus 
for  convenience,  -V  will  be  plotted  and  the  real  part  of  the 
calculated  values  will  be  compared  with  the  measured  values. 
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CHAPTER  3 


THE  METHOD  OF  SEQUENTIAL  LAYER  ADDITION 


In  this  method  of  inversion,  a  series  of  corrections  is 
made,  each  of  which  modifies  an  already-computed  curve. 

Thus,  each  correction  corresponds  to  the  result  of  interpos¬ 
ing  another  layer  above  an  original  half-space,  for  which 
V  =  0  and  whose  "exact"  resistivity  is  obtained  from  well 
logs  or  DC  resistivity  measurements.  This  method  has  the 
advantage  over  the  method  of  least-squares  curve  fitting  in 
that  one  can  observe  the  effect  on  the  curve  of  modifying  or 
adding  layers  and  thus  gain  an  intuitive  feeling  for  the 
physical  phenomena.  In  addition,  this  method  does  not  suffer 
from  the  inherent  instability  and  slow  convergence  rate  of 
least-squares  methods.  Instead  of  approximating  a  curve  by  a 
polynomial,  this  method  makes  use  of  exact  expressions. 

Rearranging  equation  (2.72), 


V  K  +  K 


n 


n 


+  V  K  + 


n 


+ 


0  (3.1) 


where 


V 

1 


V  when  n  is  even 
1  when  n  is  odd 


This  can  be  expanded  (with  sign  change,  as  explained  in  the 
last  chapter)  : 


. 


' 
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V  =  <ie-QWl  +  <2e-Q<Wl+W*> 

+ 

-Q(W!+W2+W3) 

K  3  e 

+  V<lK2e_QW2 

+ 

V<,K3e-Q(W2+W3) 

+ 

V<2<3e-QW3 

+ 

K1<2K3e-Q(Wl+W3 

It  can  be  seen  that  successive  columns  on  the  right  side 
of  (3.2)  contain  parameters  of  successively  deeper  layers; 
thus,  the  j ' th  column  does  not  contain  the  parameters  of  the 
j+l'th  layer. 

We  now  define  the  following  sequence  of  approximations: 


V1=K1e 


-QWi 


vci  =  vi 


V  -V  +  <22 


V  =V2+  <se 


-Q(W,+W2) 


-Q(W!+W2+W3) 


V. 


VC2  = 


1  +  K i K26 


-QW; 


VC3= 


V3  +  K ! K  2  e 


-Q(Wi+W3) 


1+k i k2s 


■qw2+^  <  -Q (W2+W3 ) 

-+-K  1  K  3  e 


+  k2  k  3e 


-QW. 


(3.3) 


The  VC.  can  be  seen  to  be  the  exact  expressions  for  a  (j+1)- 
D 

layered  earth  (Compare  with  the  expressions  2.75  and  2.78). 
Consider  the  residuals  defined  as  follows: 

=  vc1  -  V 

=  VC2  -  V 

=  VC  -  V  (3.4) 
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- 

FIGURE  3:  A  family  of  two-layer  curves 
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FIGURES  4(a)  and  4(b):  Three-layer  earth  models 


(a)  resistive  half-space 

(b)  conductive  half-space 
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The  process  consists  of  fitting  the  data  with  successive¬ 
ly  calculated  curves  of  VC  until  the  residual  R  is  smaller 

n 

than  the  scatter  in  the  data.  In  the  calculations,  the  real 
part  of  the  VC's  and  the  R's  is  used  and  will  be  compared 
with  the  in-phase  components  of  the  data. 

The  curve  VC^  is  a  two-layer  curve —  that  is,  the  func¬ 
tion  V  for  an  earth  model  composed  of  a  layer  of  finite  thick¬ 
ness  resting  on  a  half-space.  The  family  of  two-layer  curves 
shown  in  FIGURE  3  has  been  calculated  for  various  resistivity 
contrasts.  For  a  given  first-layer  thickness  and  resistivity, 
the  "rightmost  zero-crossing"  for  all  of  the  curves  (called 
T^)  is  the  same.  If  one  is  able  to  match  an  experimental 
curve  with  a  curve  of  the  type  VC^ ,  then  the  substructure  can 
be  interpreted  unambiguously  as  two-layered.  If  the  curve 
cannot  be  matched  with  VC-^ ,  then  a  successive  approximation 
is  calculated  based  on  the  introduction  of  a  second  layer 
above  the  half-space.  It  is  obvious  that  the  reflection  and 
transmission  coefficients  at  the  second  layer  boundary  will 
change  drastically  and  thus  from  equation  (2.70),  so  will  the 
function  V.  This  effect  can  be  seen  in  FIGURES  4(a)  and  4(b), 
in  which  resistive  and  conductive  half-spaces  respectively 
are  added  to  the  two-layer  model  for  which  k=0.5  as  shown  in 
FIGURE  3. 

In  the  process  of  interpretation,  certain  factors  must 
be  taken  into  account  which  are  part  of  the  physical  phenom¬ 
enon,  independent  of  the  method  of  interpretation.  Since  the 
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earth  is  a  resistive  medium,  there  exists  a  skin  depth,  as 
defined  in  equation  (2.25),  which  is  the  layer  thickness  at 
which  the  amplitude  has  decreased  to  1/e  of  its  value  at  the 
upper  boundary.  It  follows  from  the  form  of  the  equation 
that  for  shorter  periods,  there  is  a  correspondingly  smaller 

depth  to  which  the  em  wave  can  probe;  or  conversely,  the 

« 

effect  of  resistivity  changes  at  depth  can  have  little  effect 
on  the  shorter  period  values  of  either  V,  p,  or  Z.  Thus,  in 
calculating  the  first  approximation  to  V,  it  is  essential 
that  the  curve  match  for  the  shorter  periods  and  that  the 
departure  from  the  observed  values  be  in  such  a  direction 
that  the  next  approximation  causes  the  curve  to  match  for 
successively  longer  periods  without  departing  significantly 
from  the  good  match  already  obtained  at  the  shortest  periods. 
One  is  assisted  in  this  process  by  the  shielding  effect  al¬ 
ready  discussed.  However,  there  is  the  corresponding  problem 
that  it  requires  thicker  layers  and  higher  resistivity  con¬ 
trasts  in  order  to  be  observable  for  successively  deeper  lay¬ 
ers,  which  also  implies  longer  periods.  In  all  methods,  thin 
layers  and/or  low  resistivity  contrasts  which  lie  at  depth 
are  invisible  and  tend  to  be  averaged  with  their  neighbours. 
In  the  choice  of  parameters,  it  is  essential  to  anticipate 
the  effect  that  the  deeper  layers  will  have,  not  only  on  the 
curve  for  longer  but  also  in  the  shorter  period  region. 

FIGURES  (4a)  and  (4b)  illustrate  the  effect  of  adding  an 
additional  layer  to  the  two-layer  model  with  k  =  0.5,  for 
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FIGURE  5:  Curve  fitting  for  model  (3.5) 
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which  the  calculated  curve  is  shown  in  FIGURE  (3) .  It  should 
be  noted  that  for  the  more  resistive  substratum  (FIGURE  4a) , 
the  effect  does  not  propagate  to  the  shorter  periods  and  the 
zero-crossing  is  unaltered,  it  is  the  period  at  the  zero¬ 
crossing,  ,  which  is  used  to  estimate  the  thickness  of  the 
first  layer  (See  equations  3.6  to  3.8.).  One  then  chooses  a 
value  of  k  j  that  will  give  values  of  V  which  match  the 
straight  line  part  of  the  curve  for  values  of  T  up  to  , 
which  would  correspond  to  the  largest  value  of  T  employed  if 
the  earth  consisted  of  two  layers.  ki  determines  p2  from 
equation  (2.63) . 

In  the  case  of  a  conducting  second  layer,  as  illustrated 
in  FIGURE  (4b) ,  is  chosen  to  the  left  of  the  observed  zero 
crossing  since  the  effect  of  adding  the  conducting  third 
layer  will  cause  the  curve  not  only  to  be  pulled  down  at  the 
longer  periods  but  to  be  affected  also  at  shorter  periods, 
including  the  range  about  the  zero  crossing. 

For  the  model  shown  in  table  (3.5)  below, 


h  (km) 

R  (fi-M) 

5 

10 

30 

1000 

30 

50 

4000 

for  which  kj  =  0.818  and  k2  =  -0.635,  the  V  function  appears 
as  the  "exact"  curve  of  FIGURE  5.  If  had  been  chosen 


correctly  (h^  correct)  but 


. 


. 

■ 


V 


34 


too  large 

< 

<1 

,  then  R  would  have 

been 

0  in 

region  A. 

too  small 

> 

If 

h^ ,  h^/  and  k2  were  correct,  but 

too  large 

> 

Kl 

,  then  R  would  have 

been 

0  in 

region  C. 

too  small  A 

< 

If 

h^,  Kn  and  k2  were  correct  but 

h 2  too 

•  small, 

the  second 

approximation  curve  would  resemble  (i) ;  if  too  large,  it  would 
resemble  (ii) . 

These  considerations  apply  to  successive  approximations. 
If  the  parameters  chosen  in  the  j  '  th  approximation  (that  is, 
h j ,  k j )  are  not  correct,  then  the  (j  +  l)'th  approximation  will 
show  the  error. 

The  actual  determination  of  h^  is  made  as  follows?  The 
zero-crossing  corresponding  to  a  two-layer  earth  is  the  point 
at  which 

VC1  =  <ie"QWl  =0  (3.6) 

From  equations  (2.64)  to  (2.68),  one  can  write  the  real  part: 


k  j  e 


cos 


H 


0 


from  which 


47rh1 

/lo  PXTX 


(2n  +  1 ) t  n  —  0*  1/  2,  .  .  . 


(3.7) 


For  the  rightmost  zero  crossing,  n  =  0  and 


. 


. 
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hl  ^  10  p1T1  /8  (3,8 

which  is  identical  to  the  result  of  Cagniard  for  the  two- 
layer  apparent  resistivity  curves. 

Analysis  of  curve  shifts  due  to  this  technique  is  aided 
by  residual  analysis.  The  curve  shift  from  VC1  to  VC2  is: 


real(vc2-vc1) 


K  e~QWl  +  K, e-Q (wi+w2 ) 
REAL  [  - Z-£jLe 


1  +  Ki K2e 


-QW2 


Kle"QWM 


(3.9) 


=  REAL [ 


(1  “  k?  ) 


“Q (Wj  +  W2) 


1  +  k i k2  e 


-QW2 


(3.10) 


k  ?  ( 1  -  k 2  )ea+b  {kiK2  ebCOS (a)  +  COS (a+b) } 

1  +  2k ! k2  ebCOS (b)  +  (KlK2eb)2 


where  we  define 


a  E  -47Th1//10  pi  T 
b  E  -47Th2//10  p2  T 


(3.11) 


(3.12) 


At  T  =  T.  ,  a  =  — tt/ 2 .  For  T  <  T.  ,  a  <  -tt/2;  for  T  >  T,  , 


— tt/2  <  a  <  0.  As  T->oo ,  a->0_ • 

Defining  AR^  =  Rj+1  “  Rj /  then  from  (3.11), 

k2  (1  -  k  i  2  )  e  (b  "  TT/2)SIN(b) 


ARJ  = 

1fT=T1  1  +  2K!K2e 


(3.13) 


e  COS (b)  +  ( k i k 2 e  ) 


■/  •  '  £ 


■ 


- 
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b  =  -tt/2  t—  /Pl^p2  (3.14) 

nl 

Consider  the  case  P2  >>  Pi  (normally  h 2  ^  and  sin(b)  is 
always  negative.)  and  <2  <  0.  Ar^  >  0  and  the  effect  of  the 
correction  for  the  third  layer  pulls  the  curve  upward. 
However,  for  K2  >  0,  is  negative  and  less  significant 

than  for  k2  <  0.  For  more  complicated  models,  the  process  is 
equivalent  to  making  AR^  <  6  for  all  values  of  T,  where  6  is 
determined  by  the  scatter  in  the  data. 

In  order  to  estimate  h 2  from  T^,  one  takes  only  the 
principal  part  of  (3.10),  which  is  equivalent  to  the  process 
of  Nabetani  and  Rankin  (1969) : 


REAL  (EXP  [-Q(Wi+  W2  )  ]  ) 


Q  (Wi 


+  W2  ) 


=  tt/2 
T+T„ 


(3.15) 


For  T  <  T2 ,  (a  +  b)  <  -tt/2;  for  T  >  T2 ,  -tt/2  <  (a+b) <  0.  As 
T+co,  (a+b)-*0  .  It  can  be  seen  from  expression  (3.11)  that  if 
one  chooses  the  point  at  which  R^  =  R2 .  then  our  layer  thick¬ 
ness  would  be  given  by 


K 1  i<2ebCOS  (a)  =  -  COS  (a  +  b)  , 

or  h-  =  /10  Pg-TA{cOS~1[-iCiK2ebCOS(a)  ]  -  a} 

Ui  2  4tt 

(3.16) 

The  relationships  for  the  thicknesses  of  subsequent  layers 
would  be  more  involved.  The  numerical  approximation  involved 


'  ■ 


. 


■ 

FIGURE  6:  Choosing  f°r  model  (3.5) 


. 


« — > 


in  using  relationship  (3.15)  is  demonstrated  below.  In  the 
hypothetical  model  (3.5),  the  parameters  are: 
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«i  =  0.818 

k  2  —  ”0.635 
=  16.0  sec. 

=  41.0  sec.  ,  from  (3.15) 

If  (a  +  b)  =  ”tt/2,  then  (3.11)  becomes,  according  to  Nabetani 
and  Rankin  (1969) : 

AR  =  k1k22  (1  ^l.2)e< - ^_COS_(a) —  (3.17) 

1  +  2k i K2eDCOS  (b)  +  (K,K2eB)2 

Inserting  the  correct  parameters,  one  finds  that  at  T2  =  41.0 
seconds,  AR^  =  +0.004,  thus  verifying  that  (3.15)  is  a 
reasonable  approximation  and  as  shown  in  FIGURE  6,  causes  the 
curve  VC2  to  fit  closely  to  the  exact  curves  up  to  and  beyond 

V 

Kj  should  be  chosen  such  that  its  slope  is  equal  to  that 
of  the  straight-line  section  of  the  V  curve  in  the  region 
<  T  <  T2;  and  the  VC 2  curve  should  match  exactly  the  V 
curve  up  until  some  T  >  T2  if  the  correct  values  of  and  T2 
are  chosen.  If  not,  then 


a) 

if 

it 

does 

not 

match 

at  ,  choose  a  new  T^; 

b) 

if 

it 

does 

not 

match 

for  <  T  <  T2 ,  then  either 

(i)  T2  (that  is,  h2)  is  incorrect;  or 
(ii)  Kj  is  incorrect. 


■ 


■ 
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When  this  section  of  the  curve  is  matched,  one  goes  on 
to  the  third  approximation,  if  needed.  When  this  is  done, 
the  curve  adjustment  is: 


AR 


K  ,e-QW  4  k  ,e~Q  (W  i+W  +  k  .e'° (W  1+W  *+W  »>  .  k  ,  k  ,  e'Q  (W- +W*  > 

2  1  +  KlC!e-*  ;  KIK3e-Q(W2+W^)  7  K2K3e"QW3 


k  i  e 


-QWi 


+  k  2  e 


-Q(W!+W2) 


1  +  k  i  ic2  e 


-QW2 


K3e-Q(w,+w3+w3)  (1  _  tCi2)  (1  _  K;2) 

DENOMINATOR 


(3.18) 


In  a  similar  manner. 


Ar„  =  ^ 


-Q  (Wi  +W2  +W3  +W4  ) 


(1  -  K!2)  (1  -  K22  )  (1  - 


K32) 


DENOMINATOR 


(3.19) 


Generally,  the  curve  adjustment  between  successive  approxima¬ 
tions  is: 


ARj-l  =  [VCj  "  VCj-l] 


j-1 

k  . (e  )  n  (1  - 

^  m=l 


<  2  ) 
m 


(3.20) 


DENOMINATOR 


and  as  has  been  previously  stated,  the  process  is  continued 

until  R.  <  6,  the  scatter  in  the  data. 

3 

In  order  to  compute  the  successive  thicknesses,  one 
repeats  the  process  leading  to  equations  ( 3 . 8 )  and  (3.15): 


4tt  (• 


/10  Pi  T2  /l0  p 2  t2 


-)  =77/2 


(3.15) 


. 


J  .  .  i  !  .-A-:  r  '  rj  .  1 
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Thus , 

h2  =  /10  p2  T2/8  -  h1/p27Pi 

=  /To  T2 / 8  —  /1 0  p”i  T i  / 8  /p  2 /p  i 

=  /10  p 2/8  [  /t7  -  /t7  ]  (3.21) 

and  in  general, 

h.  =  /lO  p./8[  /FT  -  /tTT  ]  (3.22) 

1  1  1  l-l 

This  method  of  interpretation  is  implemented  using 


interactive  programming,  which  is  described  in  the  next 
chapter . 


■ 
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CHAPTER  4 

RESULTS  AND  SUGGESTIONS  FOR  FURTHER  WORK 

The  method  of  sequential  layering  was  applied  to  several 
representative  horizontally-stratified  earth  models;  the  main 
computer  program  listings  are  given  in  APPENDIX  C. 

FIGURE  7: 

As  can  be  seen,  the  third  layer  has  a  low  resistivity 
compared  with  that  of  the  second  layer.  Thus,  T^  is  chosen 
to  the  right  of  the  rightmost  zero  crossing.  The  curve  fol¬ 
lows  parallel  to  the  model  in  the  straight-line  region  of 
30-60  seconds  period.  T^  is  chosen  while  R^  is  less  than 
zero.  The  second  approximation  corrects  for  the  previous  one 
and  follows  the  curve  until  periods  of  over  1000  seconds.  At 
T^ ,  R2  is  greater  than  zero  and  the  third  approximation  cor¬ 
rects  for  the  rest  of  the  curve. 

FIGURE  8: 

As  before,  R^  is  less  than  zero  at  T^  and  .  It  should 
be  noted  that  because  of  poor  contrasts,  the  curve  has  little 
"character"  and  therefore  the  data  could  be  fitted  satisfac¬ 
torily  by  a  range  of  values  of  h ^  and  p2» 

FIGURE  9: 

The  third  layer  has  a  relatively  high  resistivity  and 
there  is  thus  an  inflexion  in  the  curve.  Since  there  is  lit¬ 
tle  resistivity  contrast  in  the  top  two  layers  and  the  second 
layer  is  thick,  T^  can  be  chosen  right  at  the  zero  crossing. 


. 


* 

. 
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Since  p3  is  greater  than  P2,  R1  is  less  than  zero  at  H  . 
FIGURES  10,  11,  and  12: 

In  all  curves  with  a  relatively  resistive  upper  layer, 
the  first  approximation  followed  the  curve  closely,  unlike 
the  previously-shown  examples.  The  previous  pattern  was 
generally  seen:  if  the  model  consisted  of  alternating  layers 
of  high  and  low  resistivities,  values  of  T  were  chosen  before 
the  residual  reached  zero;  otherwise,  after  T  was  zero. 
in  FIGURE  12  is  an  exception  due  to  the  fact  that  the  effect 
of  the  third  layer  is  small. 


It  was  found  that  the  most  important  requirement  for 
grid  matching  was  experience  with  the  graphics  display.  The 
results  shown  previously  were  far  superior  to  any  obtained  by 
other  techniques  in  use  at  the  time  of  writing.  This  method 
of  inverse  analysis  is  appealing  for  its  heuristic  merits; 
the  assumptions  of  least-squares  methods  are  unnecessary;  the 
obviation  of  complex  iterative  techniques  is  its  greatest 
actual  saving.  For  instance,  the  operator  may  match  a  fairly 
complicated  curve  quite  well  using  less  than  half  of  a  minute 
of  actual  computing  time. 

To  date,  the  techniques  for  interpretation  of  the  resis¬ 
tivity  structure  of  the  earth  have  been  based  on  the  calcula¬ 
tion  of  apparent  resistivity,  which  involves  calculation  of 
the  amplitudes  of  the  fields;  thus,  phase  information  was  not 
required  in  the  measurements.  Since  the  technique  presented 


' 


. 

. 


here  is  based  on  the  calculation  of  the  real  part  of  the 
fields,  both  phase  and  amplitude  information  must  be  known. 
In  future  field  work,  care  must  be  taken  to  obtain  accurate 
values  for  the  phase. 


FIGURE  7 :  Four-layer  demonstration  model 
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FIGURE  8: 


Three-layer  model  curve  fitting 
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FIGURE  9:  Three-layer  model  curve  fitting 
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FIGURE  10:  Three-layer  demonstration  model 
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FIGURE  11:  Three-layer  model  curve  fitting 
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FIGURE  12:  Four-layer  model  curve  fitting 
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APPENDICES 


APPENDIX  A 


Solution  of  the  Wave  Equation  for  a  Source  of  Finite 
Dimensions . 


If  the  inducing  magnetic  field  has  an  arbitrary  distribu 
tion,  the  solution  of  (2.17)  can  be  found,  following  Price 
(1962)  : 

iojt  -*■ 

E  =  e  C(z)  F(x,y)  ( A .  1 ) 

where 


F(x,y)  =  (9p/9y,  -3P/9x,  0) 


(A. 2) 


so  that  equation  (2.9)  is  satisfied  and  E  =0,  because  the 

z 

layers  of  equal  conductivity  are  parallel  to  the  surface  and 
thus  the  induced  currents  flow  parallel  to  the  surface  of  the 
earth  (that  is,  in  the  x-y  plane). 

Substituting  (A.l)  into  equation  (2.17), 


+ 


F (x,y) 


d2c  (z) 

dz  2 


=  k2elwt  ? (z)  F (x,y) 


(A. 3) 


or 


1  92F (x,y)  ,  9  2F (x,y) 

F(x,y)  9x2  9y2 


1  d2c(z) 
C  (z)  dz2 


(A. 4) 


Since  both  sides  of  equation  (A. 4)  are  independent,  they  can 
be  equated  to  the  constant  -v2.  The  resultant  equations  are 


, 


. 


+ 


+  v2F(x,y) 


0 


(A. 5) 


a2F(x,y) 

8x2 


32F(x,y) 

3y2 


and 


d2;(z) 

dz2 


(v2  +  k2)  c (z) 


(A. 6) 


where  v  is  known  as  the  "source  dimension  factor",  which  goes 
to  zero  for  plane  waves. 


APPENDIX  B 


Matrix  Computation  for  Calculation  of  Surface  Impedance 


The  matrix  M  in  equation  (2.51)  can  be  rewritten  accord¬ 
ing  to  Newman  (1962) 


b  i  ci 

ax  b2  c2 

b  3  C  3 


where 

a  i  =  1  +  k  o 
a4  =  exp[-k2z2] 


a2n-l  b2n 

c2n 

a2n 

b2n+l 

C2n+1 

a2n+l 

b2n+2 

=  exp [-k i z  i  ]  a  3 

=  ( It  K 

(B.l) 


a  =  exp[-knzn]  a2n+i=  (!+<„)  exp  [k^] 


bj  =  -  K0  b  2  =  -  Ko  b  3  =  “K  1  exp  [k  1  Z  1  ] 

b,=-KieXp[-k2Z,]...bJn+1=-Knexp[knzn]  b2n+2=-<nexp[-kn+1znl 

Cl  =  1  ,  c2  =  -1  c 3  =  (Ki  -  1) exp [-k2Zi ] 

c4=-exp [k2  z i ] . . -c2n=-exp[kn+1zn]  c2n+1= (Kn_1) exp [_kn+i Zn] 


(B .  2 ) 


■ 


. 

* 


Minors  of  this  matrix  are: 

D0  =  1 

D  i  =  bj 


Dk  =  bkDk-l  ~  ak-lCk-lDk-2  '  2  <  k  <  2n  +  2. 


(B .  3 ) 


Thus , 


D~--~  =  -K„exp[-k_J.1z_]D,_J.,  +  (1  -  k  )  exp  [  —  ( k  n  -  k_)z_]D. 


2n+2 


n 


n+1  n  2n+l 


n 


n+1  n  n  2n 


and 


(B .  4 ) 


°2n+l  =  -KnexP[knZn]D2n  +  exp[-kn(zn  -  *„_!>]  (B.5) 

for  even  and  odd  minors  respectively.  Substituting  (B.5) 
into  (B . 4 )  gives 

D0  , o  =  exp[-(k  ,,-k  )z  ]D_  -  k  exp[-(k  ,,+k  ) z  +k  z  ,  ] D0  , 

2n+2  r  n+1  n  n  2n  n  ^  n+1  n  n  n  n-1  2n-l 

(B .  6 ) 


The  determinant  det  of  the  submatrix  M^,  formed  by 
omitting  the  first  row  and  column  of  M,  can  be  treated  simi¬ 
larly  to  obtain 


D'  =  -k  exp[k  z  ] D '  , 

2n  n  e  n  n  2n-l 


exp[-k  (z 
n  n 


-  Zn-l>)D2n-2 


(B.7) 


exPf-<kn+rkn)zn]D2n  ^nexp'-(Vltkn)zn+knzn-llD2n-2 


(B.8) 


1 

~Ko 

1 

-K()exp[-k1z1]  -  K1exp[k1z1] 

exp[-(k2  -  k1)z1]  +  KQK1exp[-(k2  +  k1)z1] 


(B .  9 ) 

%  -  1 

Di  -  -Ko 

D2  =  exp[-k1z1]  +  K^exp  [k^] 

°3  =  "K0exp[“^2  “  kl)zl]  "  <1exp[-(k2  +  k1)z1] 


(B.10) 


Then  equation  (2.56)  of  the  text  can  be  written  as 


D  -  n ' 

_1  2n+2  u2n+l 

k.  D0  +  D'  . 

0  2n+2  2n+l 


(B.ll) 


* 
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PAGE  0001 


C 

c 

c 

c 


c 

c 

c 


c 

c 

c 


NLAYER  PROBLEM  IN  MAGNETOTELLUR IC S 
WRITTEN  BY  CHUCK  MGZESCN,  MARCH,  1971. 

######  h  4* *  *i  H  Hp  4ft  »  t  fe  U  tt#  ti  *  itt 

THIS  PROGRAM  CGMPUTES  AND  PLOTS  THE  VALUES  FOR  THE  TWO-,  THREE-, 
FCUR-,  AND  FIVE-LAYER  APPROXIMATIONS  AS  WELL  AS  THE  EXACT  CURVE 
FOR  AN  N- LAYER  MODEL  USING  THE  VALUES  OF  RESISTIVITIES  AND 
THICKNESSES  CALCULATED  FROM  THE  INTERACTIVE  PROGRAM  "GRAF 1C  "  . 


c 

REAL 

TT (256) , RG  AM ( 2  56 ) 

REAL 

T  (256)  ,M5)  ,  h  (  5  ) 

REAL 

REAL 

KAY  1 ,KAY2,KAY3,KAY4 

ARG ( 256 ) , ARG1 (256) ,ARG2(256) ,ARG3(256) 

REAL 

FN(256),FN 1(256), FN 2(256), FN 3(256) 

REAL  V ( 25 6 )  ,VT(256)  ,VTT<256),V7TT(25o) , VTTTT (  256) 
REAL  Vi  (  2  56)  ,V1C<  256)  ,  R1  (256)  ,R1C(  250 
REAL  V2 (256 ) , V2C( 256 ) ,R2 (256) ,R2C( 256) 

REAL  V  3 ( 2  56 )  ,V3C(256) ,R3(256) ,R3C(256) 

REAL  V4  <  2  56 )  ,V4C( 25  b) ,R4(256),R4C(256) 

CCMPLEX  V I ( 256 )  , GAM (2 56) 


COMPLEX 

CCMPLEX 

C  SORT ,C  E  X  P 

C,CK,CCMEGA,K1,K2,PI,P2,P3,P4,EX,HY 

COMPLEX 

A(256), 6(256) 

CCMPLEX 

COMPLEX 

CCMPLEX 

FNC ( 2  56 )  , V1CC (256) , FN1C (256)  ,V2CC(256) 

FN2C ( 2  56 ) , V3CC (2  56  ) 

FN3C ( 25fe ) , V4CC (256) 

DIMENSION  GAMMA <20 , DELTA (20) 


C 

C 

c 

c 

c 

c 


»»> . »>>> . >>>>>.  ..>>>>> . >»>> . »»> 


REAC(5,20C)  <  DELT  A { K ) , K  =  1,20) 

200  FORMAT  ('20  A 4.) 

READ (5,90)  NPERO 
90  FORMAT (15) 

REA  3(5,92)  (T (  I)  ,  1=1 ,NPERD) _ 

92  FORMAT ( 8F10. 2 ) 

REAC(5,93)  Tl 

R  E  A  0 ( 5 , 9  3  )  KAY1 _ 

E  A  C  (  5 , 9  3  )  T  2 
READ(5,93)KAY2 


_ REAC( 5,93  )  T  3 

RE  AC (5, 93) KAY  3 
R  EAD ( 5, 93  )T4 
REAOtS  »9  UKAY4 
3  FORMAT (F10 .2 ) 

C 

C _ _ _ 

•  R  E  A  C ( 5 , 9  0 )  NLAYER 


C 


RE  AD ( 5 , 91 ) ( H(  I  ) , R (  I ) »  I  =  1, NLAYER) 


. 


i  I  ,  •  )  *  • 


) 
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G  COMPILER 


C 9-17-71 


10:45.02 


C 

c 


91  FORMAT ( 2  F 10 . 1 ) 


C 

C 

C 

C 

c 

c 


FIRST,  WE  CALCULATE  THE  VALUES  OF  "V"  .  .  . 

IN  THIS  SECTION,  WE  CALCULATE  THE  "V"-CURVE.  THIS  PROGRAM 
SECTION  IS  BASED  ON  THE  CONTINUITY  OF  THE  ELECTRIC  AND 
MAGNETIC  FIELC  COMPONENTS  AT  THE  BOUNDARIES  BETWEEN  THE 
HORIZONTAL  LAYERS  OF  THE  EARTH. 


C 

C 

C 

C 

c 

c 


C  IS  THE  SQUARE  ROOT  OF  -1« 

TP  I  IS  2*PI 

A(NLAYER)  =  0.  AND  B (NLAYER)  =  i.,  CONSISTENT  WITH  THE  FACT 
THAT  IN  THE  HALF-SPACE,  THERE  CAN  BE  NO  REFLECTED  WAVE  AND 
THE  TRAVELING  WAVE  MUST  DECAY.  THE  BULK  OF  THESE 
CALCULATIONS  CCULD  BE  USED  FOR  CALCULATING  THE  APPARENT 


C 

c 

JL 

c 

c 

c 


RESISTIVITY  AND  PHASE  ANGLE  OF  THE  cLECTRIC  FIELD  WITH 
RESPECT  TC  THE  MAGNETIC  FIELD.  AFTER  STATMtNT  3,  WE  HAVE 
INSERTED  THE  CALCULATION  FOR  "V" ,  THE  MODEL  CURVE. 

THE  PARAMETERS  "  ICU"  AND  "NUPLQT*  WHICH  APPEAR  IN  ALL  OF  THE 
MOZESCN  NLAYER  PROGRAMS  ARE  PLOTTING  PARAMETERS  USED  IN  THE 
SUBRQLTINE  KNOWN  "GPL"  AND  WRITTEN  BY  R.  F.  HOWELLS 


C 

C 

C 

c 

c 

c 


OF  THE  CIVIL  ENGINEERING  DEPARTMENT. 

THE  ALPHANUMERIC  DATA  READ  IN  UNDER  THE  NAMES  "ALPHA*1, 
"BETA'S  "GAMMA'S  "PELT A"  ARC  USED  FOR  TITLING  THE  PLOTS  DONE 
ON  THE  CALCOMP  PLOTTER. 


c  =  (0.0, 1.0  ) 

TP  I  =  6.  2831853 


CG  2  IJ  *  1 , N  PERD 
COMEGA  =  C*TP  I/T  (  I  J  ) 
CK  =  -0.2*TPI*CGMEGA 


TT(  I J ) =SQRT ( T  (  I J  )  ) 

A (NLAYER )  =  (  C  .C, C.O  ) 
B( NLAYER)  =  (1.0, 0.0) 


NL1  =  NLAYER  -  1 


DO  3  J  =  1 ,  N  L 1 
I  =  NL1-J  +  1 
HI  -  H(  I  ) 

K 1  =  GSQRT (CK/R(  I  )  ) 

K  2  =  CSQRT ( C  K  /R (  1  +  1 )  ) 
PI  =  (Ki  +  K2)/( 2.*K1  ) 


P2  =  (K1-K2)/(2.*K1) 
P  3  =  -K1*H1 
P  4  =  K1*H1 


A  (  I  )  =  P1*A { 1+1 ) *CEXP( P3)  +  P2*8(I  +  l)#CtXP(P3) 
B  ( I  )  =  P2*A  (  I  +1  )  *CEXP(  P4  )  +  Pl*a(I+l)*CtXP(P4) 


3  CONTINUE 


=  REAL ( A( 1  )/B ( 1  )  ) 


V  (  I  J  ) 


r  * 


) 


l  V 

2 '  |j 


■ 


4  *  1 


f  *  l* 

*  o  )\(  .  *-.  )  =  ^ 


-  1  I 
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C 


C 

C 


C 


2  CONTINUE 

WRITE (6,201)  NLAYER 

201  FORMAT ( 40X  ,  1 2  ,  *  LAYFRED  ISOTROPIC  EARTH*/) 

WRITE! 6, 202  )  ( M  I  ),R( I  ), I  =  1, NLAYER) 

2  02  FORMAT ( 4  0 X , F 1 0  «  2 ,  5  X  ,  F  1  j  .  2  ) 

WRITE(o,205) 

205  FORMAT!//) 

WRITE! 6,204 ) 

204  FORMAT <2 <2X, • T(SEC)  1  ,9X,»SQRT!T) *  ,8X,  ’GAMMA ' ,  12X 
WR I TE (6 , 205 ) 


WRITE!6,203) !T(M) ,TT(M) ,  RG AM( M) ,V!M) ,  M=1 , NPERD  ) 
203  FORMAT { 2 ( F10 . 3 , 5X , FI C .3 , 5X , F 10. 3 , 5X , F 10.  3 ) ) 


IOU  =  8 
MJPLOT  =  1 

CALL  GPL! T  ,V ,V, NPERD, NUPLUT , 2 , 1 , 5, 3 ,0. , ^ , 9  * , 

C- • 2 , • 24 , 5. , DELTA, ICU) 

10  CONTINUE 


•  V  *  , 5X)  ) 


C 

c 

c 

c 

c 

c 


sjc  si;  sjs  3js  s!;  sj:  #  sjc  sjc  sjs  sjc  s[:  sjc s;c  sjc  sjc  sjc  sjc  sj;  s^;  sjc  sj;  sjc  ;J;  sjc  sj;  sjc 
sjc  s;c  sjc  sjc  s}c  sjc  sjc  sjc  sjs  sj:  sjc  sjc  sj t  s|c  sj;  s£  sjc  sjc  sjc  #  sjc  sj;  s';  sj;  s^  sj;  sj; 


sO  UL  O/  x  O-  %•/  *■  <  V-  vf  %■<  «i«  «V  •<#»  ■*>«•  v1.  y.  y.  »'»  >'<  v#  s'#  \*«  yu  y»  v-  v-  x 

^  /,*  ^  ^  «y*  •'j*-  r,k  *|»  #,  ■'(S  <Y"  *,»  “Y"  'r  •'(»  -r*  6*  4'  '■*  '»■*  ^  4*  t 

sjc  sjc  s,1;^;  *  sjc  s';  sj:  sjc  sjc  sj;  sj;  sjs  sjs  sj:  sjc  s^  s£  sjc  sjc  sj;  sj;  sj;  sj;  sjs  sjc  sj;  sj;  sj; 


FPI  =  12*5663706 
S IG  =  1 ./SORT! R( 1  )  ) 

HI  =  SQRT(R( 1)*T1*10. )/8. 

_ HH  =  FPI*H1*SIG _ _ _ 

C 

c 

WRITE! 6, 1371 ) hi 

13  71  FORMAT!  10X, ■ HI  =  ',F8.3//) 

C 

C _ _ _  _ 

C  WE  NOW  9EGIN  CALCULATION  CF  THE  TWO-LAYER  APPROXIMATION, 

C  WHICH  DEPENDS  ON  THE  C6SERVED  VALUE  OF  Tl,  THE  "RIGHT- 

C  MOST  ZERG-CRJSS INC  OF  THE  TIME  AXIS  BY  THE  FUNCTION  V. 

C  MOW  THAT  wE  HAVE  THE  VALUES  FOR  THE  MODEL*  WE  SHALL 

C  CALCULATE  THE  VALUES  OF  THE  FIRST  APPROX  I  MAT  I  ON , 

C _ WHICH  IS  PS SENT  I  ALLY  A  2-LAYER  CURVE.  THIS  CURVE _ 

C  SHOULD  FIT  THE  MODEL  CURVE  UP  TO  A  PERIOD  T2. 

C  THE  CALCULATIONS  ARE  STRAIGHTFORWARD  AND  EXACTLY  AS 

C  LAID  OUT  IN  T  H L’  MCZESON  THESIS* 

C 

C 

c _ _ _  - 

•  DO  1300  J  =  i , NPERD 
ARG ! J )  =  -HH/SQRT !10.*T( J ) ) 

F  N ( J )  =  EXP! ARG! J)  )*CCS( ARG!J ) ) 


t  )  I 

•  T  t  '  l 


♦  £  .hi 


M  (  4  >  .  .. 


. 

. 


■  I. 


ooooooooo  no  o  ,  o  noon 
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FNC(J)  =  EXP ( ARG { J )  ) *  *C*S I N ( ARG { J )  ) 
V1CC(J)  =  KAY  1* ( FNC  <  J  )  +  FN(J)) 

VI ( J  )  =  KAY1*FN( J ) 

_ V  1  (  J  )  =  KAY1*FN(  J) _ _ 

V1C(J)  =  VI ( J ) 

RKJ)  =  V  1  (  J  )  -  V  (  J  ) 

R1C  {  J  )  »  VIC ( J)  -  V ( J ) 

1300  CONTINUE 

WRITE (6»1362)T1 
1362  FORMAT!  IPX,  «  T1  =»  ,F6.2//) 


1361  WRITE<6,1360 )  KAYi 

1360  FORMAT  (  IX  ,  1  VALUES  FOR  THE  FIRST  APPROX  I  MAT  ION*  WITH  KAYi  =  MX, 

CF5.2/) 

WRITE(6,2C6) 

206  FORMAT ( 22X , * T< SEC » , 9X, • SORT (  I )  •  ,  10X  , » VI • , 8X, • R1 «  *  12 X, • VIC » , 12X , 
C 'RIC*  ) 


WRITE(6,205) 

WRITE(6,213)<T(M),TT{M),V1<M),R1(M),V1C(M),R1C{M),M=1,NPERD) 
213  FORMAT ( 20 X, F 10*  3 , 5X , F1Q. 3 , 5X , FI  0, 3 , 5X, F 10,3 1 5X , F10«3 , 5X> F 10.3 > 


WRITE{6,205) 


NUPLOT  =  2 

CALL  GPL ( T , VIC  ?V1C  ,  NP  ERD , NUPLOT ,2,1 , 5,3 ,3.  ,4.  ,9., 
Q»2*«24>5«,  DELTA,  IGUJ 
NUPLOT  a  3 

CALL  GPL  I T,R1C ,V,NFERD, NUPLOT ,2, 1, 5,3,0. ,4. ,9., 
C-.2  ,.24,5.  , DELTA,  IOU) 


NEXT,  WE  CALCULATE  THE  THREE-LAYER  APPROXIMATION,  WHICH 
IS  DEPENDENT  UFCN  T«E  POINT  AT.  WHICH  TFb  FIRST  APPROX- 
I  MAT  I  ON  LEAVES  THE  V-CURVE,  T2. 


FRAC  =  ( 1  .-KAYI )/{ 1.+KAY1 ) 

R  (  2  )  =  R  (  1  )  *{  (  KAYi  +  1 .  )  /  (  1.  -  K  AY  1 )  )**2 
_ H2  =  1. /FRAC* ( SORT ( 10«*R(  1  )*T2)/8.  ~H1 ) 

C 

WRITE (6,990 )R( 1) 

990  FORMAT  (  IX,  '  R(  1  )  =  SF6.1//J 
C 

WRITE( 6, 1410  )  H2 

1410  FORMAT  ( IX,  1  H2  =  SFb.2//) _ 

•  WRITE(6,i400)R(2) 

1400  FORMAT ( IX , ’R ( 2 )  =  *,F6.1//) 

C 


i  (  J>  )  T  I 


«  ) 


*  f  '  ’  f  '  «  »  *  :  '  f  T  '  (  i  )  I  »  ,  *  '  >  ,  »  t  .  >  T  iOn 


»  1 )  I  -  1 


. 


. 


If)  »  f  > 


t  *  ■  i  )  1  , 
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WRI T  E  (  6  , 1363  )T2 
1363  FORMAT ( 5X , *  T 2  =',F8.2//) 

,  DO  1330  J  =  1 , NPERD 

X. - ARG1(J)  =  Aftfil.l  )*FRAC*H7/m _ _ _ 

FNKJ)  =  EXP(  ARG1  {  J  )  >  *CGS  (  ARG  1(  J  )  } 

FNIC(J)  =  EXP!  ARGH  J  )  )*C*SIN(  ARG1IJ  )  ) 

VTT ( J)  =  KAY2*FN1 (J) 

V2(J)  =  V1C(J)  +  VT T ( J ) *Fh( J ) 

V2CC !  J  )  =  (V1CC ( J ) +KAY2*!FN1 ( J )+FNlC(  J)  ) *(  FM  J  ) +  FNC ( J ) )  )  / (  l.  +  KAYi* 
I _ C  K  AY  2*  ( FN 1  (  J  )  +  Ffsi  1  C  (  J  )  )  I _ 

V 2C ( J )  =  REAL { V2CC ( J  >  ) 

R2U)  =  V2(  J)  -  V  <  J  > 

R2C (J)  =  V2C ( J )  -  v ( J ) 

1330  CONTINUE 
C 

I _ WRI TE! 6,205) _ 

C 

1341  WRI T  E  (  6  , i 340 )  KAY2 

1340  FORMAT ( i X , *  KA Y2  ;'iX,F5.2/) 

C 


WRITE(6,1370 ) 

1370  FORMAT ! IX , 1  THE  VALUES  FCR  THE  SECOND  APPROXIMATION  ARE : *  //  ) _ 

WRITEI6, 1321 ) 

1321  FORMAT ( 2 5X , *  T ( SEC ) ‘ , 9X , *  SORT ( T )  »  ,  10 X  ,  ' V 2  * , 1 5X , ' R2 1 , 12 X,  ' V2C • , i2X , 
C ' R2C *  ) 

C 

WRI T E (6  f  205 ) 

_ WRITE(6,1322)!T(m),TT(M),V2(M),R2(M)  ,  V2C <  M) , R2C( M ) ,M=1 , NPERD) 

1322  FORM AT ( 20 X, F 10. 3, 5X, FI 0.3, 5X, FI  0.3 , 5X, F10.3, 5X , F10. 3 , 5X, F10.3 ) 

WRI TE ( 6,205 ) 

C 


C 


NUPLOT  =  4 

CALL  GPL<  T, V2C, V  ,NPERD#  NUPLGT ,2, 1, 5,3,0. ,4. ,9.  , 
C-.2,.24,5.,DELTA,ICU) _ 

NUPLOT  =  5 

CALL  GPU  T  ,R2C  ,V  ,  NPERD ,  NU  PLOT  ,  2 , 1 , 5, 3 , 0  .  ,  4  .  ,  9  .  , 
0.2, .24*5. .DELTA . IOU) 


$$$$$$$$$$$$$$$$$$$$ 

*  :]£  #  ajc  %  #  ^  1?  #  £  *  *  *  >}:  &  #  *  ❖  *  *  #  *  *  #  *  V  ?i<  »  »  »  »  »  ft  ft  ft 


$$$$$$$$%$$$$$$$$$$$ 
^  O/  V/  v1/  «,</  nI.  «,v  »*<  v*  o/  % 

'j'  'I'  /p  #(“•  rp  zp  /p  /p  /|S  /p  »p  ^p  *,v  -'I'*  /p  /j\ 


NOW,  WE  CALCULATE  THE  FOUR-LAYER  APPROXIMATION,  USING  THE 
OBSERVED  VALUE  OF  T3. 


FRAC1  =  11.  -  K A Y 2 ) / (  1 •  +  KAY2) 

R { 3 )  =  R!2)/(FRAC1)**2. 

H3  =  SQRT(R(3)/R(1)  >  *  (  S_QRT  <  10  .  *R  <  1  )  *T3  >  /S.  -  FPAC*H2  -  Hi) 


WRITE(6,1530)F3 

1530  FORMAT ( IX, * H3  =  *  ,F6.2//) 

•  WRI TE (6, 152  J  )P ( 3  ) 

1520  FORMAT! IX, *R!3)  =  *,F6.1//) 
WRITE(o,  1540  )KAY3 


v  )  LJLHM _ i 


»  •  b .  y.  *,xi )  t/.  a  »  u*n 

t • i  •  o 

»  )  r  t 

t  . 


»  .  *  1  *  »  ,  «  t  ,  }  * 


. 


U(l  1  ,c);  T  t$l4 

- 


ooo  o  o  o  o  ooooooooooo 
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1540  FORMAT ! IX , 1  KAY < 3 )  =  * ,1X»F5.2/) 
WRITE! 6 , 1364  )T3 
1364  FORMAT  < IX  » ' T  3  =  SF8.2//) 

£ 


10 : H6 .02 
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_ _ DO  1510  J  s  ItNPERO 

ARG2 ( J  )  =  ARG1 ( J  )  *FRAC 1+H3/H2 
FN2(J)  =  FXP ( ARG2 ( J ) ) *CGS { ARG2 ( J ) ) 

_ FN2C  !  J  )  =  EXP(ARG2( J  )  ) *C»SIN! ARG2! J )  ) _ 

VTTT(J)  =  KA Y3Y  FN2 ! J ) 

V  3  ( J  )  =  V2C { J  )  +  VTTT!  J)*FN(JI*FNH  J) 

V3C  C  (  J  )  =  ( V1CC  (  J)  +  RAY  2* ( FN1  !  J  )  +  F  N 1C ( J  )  )  *  !  F.N  (  J  )  +  FNC  <  J  )  )  + 
CKAYJ*  (  FN1  (  J  )  +  FMC(J))=MFN(J)  +  FN  C  (  J  )  )  *  (  FN2  (  J  )  +  FN2C(J)) 

C  +  V1CC { J ) *KAY2*KAY3* (FN2( J  )  +  FN2C ( J ) ) ) / ( 1 . 0  +  K AY  1 *K AY2 * 

_ C!_FN1!J)  +  FN1C(J))  +  KAY 1*K AY3* ( FN 1 ( J )  +  F  NIC  (  J  )  )  *  (  FN2  ( J  )  +  FM2C  (  J  )  ) 

C+KAY2*KAY3* (FN2< J)  +  FN2C(J))) 

V  3C  !  J  )  =  REAL (V3CC( J )  ) 

R3  !  J  )  =  V  3 ( J  )  -  V  !  J  ) 

R 3C { J )  =  V3CC  J)  -  V <  J 1 
1510  CONTINUE 


WRITE(6* 1540 )KAY3 


WRITE!6,1550) 

1550  FORMAT ( IX ,» THE  VALUES  FOR  THE  THIRD  APPROXIMATION  ARE:'//) 


WRITE! 6,1505  ) 

1509  FORMAT! 2 5X,  *  T( SEC  )  1  ,9X, *  SQRT( T)  * ,  10 X,  *V3» , i 5X , 'R3'  , 12 X, • V3C •  ,12X, 
0^30 


WRITE! 6,205) 

W  R I  T  E  !  o  ,  I  50  8  )  !  T  !  M  )  ,  T  T  (  M  )  ,  V3  (  M  ) ,  R3  !  M  ,  V3  C  (  M  )  ,  R  3C  (  M  )  ,  M-  1 ,  i\  P  E RD  ) 
1508  FORMAT !20X, F i C . 3 , 5X , F 1 0. 3 , 5X , F 1 0 . 3 , 5X , F 1 0. 3 , 5 X , F i 0 . 3 , 5X,F10.3) 


NUPLCT  •=  6 

CALL  GPL(T , V  3C , V3C »  NP ERO, N LPLOT  ,2,1  ,5,3,3. ,4. ,9. , 
C-.2,.  24,5.  , DELTA, IOU) _ 

NUPLOT  =  7 

CALL  GPL (T ,R3C,R3C,NPERD, NUPLOT, 2, it  5,3, 3. ,4. ,9. , 
C  -  •  2  ,  .  2  4 , 5  •  ,  D . :  L  T  /\ ,  T  C  U  ) 


»>>>, 


»»> 


>>>>> 


>>>>>. 


>>>>>, 


>>>>> . »>>> 


FINALLY,  WE  ShALL  CALCULATE  THE  FIVL-LAYER  APPROXIMATION, 
USING  THE  OBSERVED  VALUE  CF  T4 ,  THE  POINT  AT  WHICH 
THE  THIRD  APPROXIMATION  LEFT  THE  V-CURVE. 


FRAC2 


( 1.-KAY3  )  /(  1.+KAY3) 


* t  >i  nm 


)  (  (U  J  ♦  I!  \  )  .) 

t  ) 


i  ,ci)3T  1 54W 


*  i  t(  )  7 1  *w 


* 


»  «  .  *  I  f  I  , 

» 

•  t  <■  “ 


» •  ' «  .  * 


;>  '  %  -  tU  . .  ,  J  |a„U  1 

. 


■ 
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R  ( 4  )  =  R ( 3 ) / { FRAC  2  )  **2 . 

H4  =  SQRT(R ( 4 )/R ( 1 ) )*( SQRT ( 10. *R ( 1 ) *T4I /8.  - 
C-FRAC*FRAC1*H3  -  hi) 


FR AC*N2 


nn  1610  J  =  1 ,  NP  ERD 

ARG3IJ)  =  AR  G2  (  J  ) #F  RAC  2*H4/H  3 

FN3 (J)  =  EXP  ( ARG3( j ) ) *CGS( ARG3( J ) ) 

FN3C ( J  )  =  EXP( ARG3( J ) )*C*SIN{ ARG3( J ) ) 

VTTTT(J)  =  KAV4»FN3!J) _ _ _ 

V4(J)  =  V3C  (  J  )  +  VTTTT  (  J  )  *FN<  J)  *FN1  ( J  )*FN2(  J  ) 

V4CC ( J )  =  (VICC(J)  +  KAY2 /K  AY  1* V  ICC ( J ) *( FN1 ( J  )  +  FN1C(J))  + 
CKAY3/KAY1*V1CC(J)*(FN1(  J)  +  F  N 1 C  (  J  )  )  *  (  F  N  2  (  J  )  +  FN2C(J))  + 
CKAY4/KAY 1*V1CC ( J ) * ( FN1 ( J  )  +  FN1 C < J ) ) * ( FN2 ( J )  +  F  N  2  C  (  J  )  )  *  ( F  N  3  (  J  )  + 
CFN3C ( J ) )  +  KAY2*KAY3*V1CC(J)*(FN2<J)  +  FN2C(J))  +  V1CC(J)* 
C(FN2(J)  +  FN2C(  J  )  )*(  FN3(  J  ?  +  FN3C  (  J  )  )  *KAY2*KA  Y4  +  ViCC ( J ) *KAV3* 

CKAY4* (FN3 ( J)  +  FN3C(J))  +  KAY2*KAY3*KAY4/KAY 1*V1CC ( J ) * (FN1C( J ) 


C  +  FNKJ)  )*<FN3(J)  + 
C+KAY1*KAY3* ( FN1  ( J  > 
C (FN2( J )  +  F  N 2  C  (  J  )  ) 
CFN2C! J) )* (FN3( J)  + 
C  ( F N 3  (  J  )  4-  FiBC(J)  ) 


FN3C  (  J  )  )  )  /  (  1  •  +  KAY2'‘£KAV1*(FN1(J  )  +  FMC(J)) 
+  FNiCU  )  )*<FN2(  J)  +  FN2C  (  J  )  )  +  KAYc*KAY3v 
+  KAY1*KAY4*(FN1(J)  +  FiMl  C  (  J  )  )  *(  FN2  ( J  )  + 

FN3C ( J ) )  +  KAY2*KAY4* (FN2( J)  +  FN2CIJ))* 

+  KAY3*KAY4» (FNS ( J )  +  FN3C(J))  +  KAYi»KAY2* 


CKAY3*KAY4*( FN1( J )  +  FN1C ( J ) ) * < FN 3 ( J ) 
V4C< J)  =  REAL ( V4CC { J ) ) 
y  (  j )  -  V4(  J  )  -  V  (J) 

R  4  C  (  J )  =  V4C  (  J )  -  V(J) 

1610  CONTINUE 


+  F  N  3  C  (  J  )  )  ) 


C 


1368 


1620 


WRITE (6, 1368 )T4 

FORMAT (10X, ' T4  =*,F8.2//) 


WRITE(6 , 1620 )R(4 i 
FORMAT  (  L  X  .  *  K  (  4  )  =  SFb.l//) 


WRITE (6, 1630  )F4 
1630  FORMAT!  IX,  1  H4  =  '  »F6.2_//  ) 


C 

C 


P { 5 )  =  R !4)*(i.  +  KAY4  )  **2  .  /  ( 1 »~KAY4 ) **2« 


1791 


WRITE (6, 1791 )R< 5 ) 

FORMAT { 20X , ' R ( 5 )  =*,F10.3//) 


RES  =  0. 

DO  8888  J  =  1  ,  NP  ERD 
RES  =  RE  S  +  (R4C< J )  ) *»2 


— 


8888  CONTINUE 

SO  =  SQRT(RES/( NPERD-1.  ) ) 

WPITF(u,h889 LSD 

•  FORMAT (20Xf • STAND ARO  DEVIATION  -*»F10.5//I 


C 

C 


■  WRI  TE-  (6  ,  1640  )KAY4 

1640  FORM AT ( 1 X VALUE S  FOR  THE  FOURTH  APPROXIMATION,  WITH  KAY4-  *,F5.2/ 
C/)  _ 


)  X  1  (  h  )  t\t  '  >  >  Tyf  ; 

I  <  -  l- 


(  (  L  \  )< 


U  I  !  L11IVL 


'  \  . 

)  .  )  - 


(1  )+* 


) 


tv)  (  1-v  Ji  tc  >  -}  T  I  1 


n  o 
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C 

WRI TE (6,205 ) 

WRITE! 6,1367) 

13  67  FORMAT (25X, *  1  T  (  SEC  )  1  ,  9X ,  1  SORT  <T)i,UX<,V4,<15Xt  1  R  4 »  » 1  ^  X  *  '  V4C  •  ,  12X. 
C  *  R4C  *  ) 

WRITE! 6 ,205) 

WRI TE (6 , 1508 ) ( T{ M ) , TT ( M) , V4! M) , R4( M ) , V4C < M ) , R4C ( M ) , M= 1 , NPERD ) 


I  _ NUPLOT  =  8 _ _ 

CALL  GPL ( T  ,V4C,V4C,NPER0, NUPLOT ,2,1, 5, 3, 3,, 4. ,9., 
C-. 2  ,. 24,5* , DELTA , ICU ) 

_ _ _ _  NUPLOT  f  -9 

CALL  GPL (T ,R4C,R4C, NPERD, NUPLOT, 2, i, 5,3, 3. ,4. ,9., 
C  —  «  2  ,  •  2  4  ,  ,  0  n  L  T  A  ,  I  l  U  ) 

I _ STOP _ _ _ _ _ _ 

END 

MEMORY  REQUIREMENTS  01CFEA  BYTES 
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C 

c  ##########* *###  H#  nit  II  tint*#  HUH##*  It  MW##*###  §*#**#  ############*##  4  It 

C  NL AYER  PROBLEM  IN  MAGNE TOTELLUR I CS 

_ £ _ WRITTEN  BY  CHUCK  MCZESCN,  MARCH  »  1971. _ 

C  ADAPTED  FOR  INTERACTIVE  USE  BY  STEVEN  HOSKIN,  JUNE  1971. 

C  *  ######  ###  MHHt  44444444444444444444444444 

C  THIS  INTERACTIVE  PROGRAM  ASSISTS  THE  GRID  USER  TG  MATCH  TWO-.* 

C  THREE- ,  FOUR-  AND  FIVE-LAYER  APPROXIMATION  CURVES  TO  THE  EXACT 

C  CURVE  FOR  AN  N-LAYER  MODEL. 

_ c _ 4444444444444444  4444n44 4444 444444^44444444^444 _ 

c  44 44 444444 H44444444444 44444444444 44444444 4 H44 4444444 

C 

_ _  SUBROUTI  J E  GR AF IC 

LOGICAL*!  ERR/. FALSE. /, TR/.TRUE* /, F/. FALSE./ , GN/, TRUE. / 

INTEGER*2  AX ( ICC > , A Y ( 100 ) , STRING ( 10 ) 

_ INTEGER  BLK(  8)  ,  STATUS,  KEY _ _ _ _ _ 

REAL  RGAM< IOC  ) , T< 100 ) , R (10) ,H( 10) 

REAL  KAY1 ,KAY2 ,KAY3 ,KAY4, KAY5 

REAL  ARG{  LOO ) ,ARG1 ( 100) , ARG2( 10Q )  ,ARG3( 100  ) 

REAL  FNC100 ) , FN1 ( 100 ) , FN2( 100) , FN3( .00) 

REAL  V( 100) , VTUOO) ,VTT(100) ,VTTT{ 100 ) , VTTTT (100) 

_ REAL  VI (  100)  tViCC 100) ,V2( 103) , V2C< 100) , V3 ( li 0 ) , V5C ( 1 00 ) _ 

REAL  V4 (100  )  ,  V4C( 100) ,R1( 100 ) , R1C( 100) ,R2« 100) ,R2C( 100) 

REAL  R3<  100  )  ,R3C ( ICO) , R4  (100) ,R4C( 100) 

COMPLEX  VI  (100)  ,GAM( 100) 

COMPLEX  CSQRT  ,CEXP 

COMPLEX  C,CK ,C0MEGA,K1, K2, PI, P2, P3,P4,A( 100)  ,B(  100  ) 

_ CCMPLEX  PNC  ( ICC  )  ,  V1CC  (100  ) ,  FiMIC  (  100  )  ,V2CC(  100) _ 

COMPLEX  FN2C (  iOO )  ,  V3CC ( 1 00  )  , FN3C ( 100 ) , V4CC( 100 ) 

C  >>>>> . »»> . >>»> . >>>>> . >>>>>.  .  . .  .  >>»> . »>>> 

C 

c 

c _ 

C  SUBROUTINE  SCREEN  INITIALIZES  THE  GRID  DISPLAY 

C 

CALL  SCREEN 
C 
C 

C _ NPERD  IS  THE  NUMBER  OF  PC  I  NTS  FOR  WHICH  DATA  WILL  BE  CALCULATED 

C 

10  WRITE(6,80) 

_ 80  FORMAT ( *  ENTEF  NPJEJjt D  •  ) 

READ (5,90)  NPERO 
90  FORMAT (15) 

C 

c 

C  FILE  1  CONTAINS  TEE  TIMES  FOR  WHICH  VALUES  OF 

C  V  WILL  BE  CALCULATED 

C 

RE AD ( 1 ,92 )  ( T(  I )  ,  I =1 , NPERC) 

_ 92  FORM AT ( 8  F  10 . 2  ) _ 

•  GO  TO  16 

15  CALL  DLPEN(  2,  I  X,  I  V,  I  TYPE  ,  ID,  BLK  ,E16) 

GO  TO  17 


' 


t  f 

* 


•  ;  '  '  *  ■  •  <  <  •  •  •' . ■  (<•:**..  .<<<  < 


* 


■ 


. 
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C 

C 

C 

£ 


C 

c 

c 

c 


Nl. AYE  R  IS  THE  NUMBER  OF  LAYERS  IN  THE  EARTH  MODEL 


16  WRITE(6,33) 

83  FORMAT ( 1  ENTER  NL AYER ' ) 
READ! 5,90)  NLA YER 


IPGSN  INDICATES  HOW  SCRtEN  COORDINATES  MUST  BE  ESTABLISHED 


WRITEI6,  111) 

111  FORMAT  (  *  ENTER  1  IF  RESISTIVE  CURVE;  -1  IF  CONDUCTIVE  CURVt1) 
READ  (5,90)  IPGSN 
IFUPCSN  .EQ.  1)  I  HER  E  =  100 
IF ( 1  PCS N  .EG.  -1)  I  HERE  =  924 


ESTABLISH  RANGE  ON  X-AXI.S 


WRITE(6,112) 

112  FORMAT ( 1  ENTER  LCwFST  VALUE  ON  X-AXIS* ) 


RE  AD ( 5 , 92 )  XSCALE 
DO  50  I  =  1 1  N P ERD 
50  T(I)  =  T ( I )  *  XSCALE 


IN  THIS  SECTION  WE  CALCULATE  OR  READ  IN  VALUES  OF  "V”  .  .  . 

THE  CALCULATION  CF  THE  "V^-CURVE  IS  BASED  ON  THE  CONTINUITY  UF 
THE  ELECTRIC  AND  MAGNETIC  FIELD  COMPONENTS  AT  THE  BOUNDARIES 
BETWEEN  THE  HORIZONTAL  LAYERS  OF  THE  EARTH. 

C  IS  THE  SCUARE  ROCT  OF  -1, _ 

TPI  IS  2* P I 

A(NLAYER)  =  0.  AND  B  ( NL  AYER. )  =  1.,  CONSISTENT  WITH  THE  FACT 
THAT  IN  THE  HALF-SPACE,  THERE  CAN  BE  NO  REFLFCTED  WAVE  AND 
THE  TRAVELING  WAVE  MUST  DECAY.  THE  BULK  OF  THESE 
CALCULATIONS  COULD  BE  USED  FOR  CALCULATING  THE  APPARENT  ' 
RESISTIVITY  AND  PHASE  ANGLE  OF  THE  tLECTRIC  FIELC  x I T H _ 

RESPFCT  TO  THE- MAGNETIC  FIELD.  AFTER  STATMENT  3,  WE  HAVE 
INSERTED  THE  CALCULATION  FOR  "V",  THE  MODEL  CURVE. 


IDAT  =  0 


■  C _ _ _ - _ 

C  INDICATE  IF  DATA  CURVE  TO  BE  SIMULATED  OR  REAL  DATA 

C 

■  WRITE! 6,113) 

113  FORMAT! *  If  V  IS  TO  BE  SIMULATED,  ENTER  i') 

R E AD ( 5 , 90  )  IDAT 

_L _ IF  (IDAT  .ME.  1)  GO  TC  250  _ 

C 

c 

C  H  AND  R  FOR  EACH  LAYER  IS  ENTERED 


S  T  »>TVMI 


I  ,  I  , 


' 


* 


. 


. 
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DO  100  I  =  It  MAYER 
WRITE (6t84)  I 

84  FORMA  T (  1 _ ENTER  H  AGO  R  FOR  IAY~r  «  .  j  i  ) 

100  RE*D(5,91 )  H(I),R(I) 

91  FORMAT (2  FI 0. 1 ) 


C 

250  C=  <0.0, IcO) 

1 _ TP  I  =  6.2931S55 _ 

FPI  =  2 • 0*T  P  I 
C 

DO  2  IJ  =  1,NPERD 
COMEG A  =  C*TPI/T(IJ) 

CK  =  -0«2*TPI*CCKE.GA 
I _ IF  (I  DAT  .EQ.  0  )  GO  TO  2 

A(NLAYER  )  =  {0oC,0.0) 

8 ( NL  AY  ER )  =  <  1. 0,0.0) 


NL1  =  NL AYER  -  1 
CO  3  J  =  1 , NL1 
I=NL1-JU _ 

H 1  —  H  (  I  } 

K 1  =  CSQRT ( CK/R  <  I  )  ) 

K2  -  CSQRT ( CK/R {  1  +  1)  ) 


R E AD ( 5,91  )  R (  1  ) 


C  INSERT  FACILITY  TO  READ  DATA  HERE 

85  FORMATt  10F10.3  ) _ _ _ 

200  SIG=1.0/SQRT(R( 1) ) 

WRI TE ( 6, 1 01 ) 

101  F  0  RM AT <  2  1H  DATA  CURVE  DISPLAYED) 

C 

C  CISPLAY  DATA  CURVE 

C _ _ _ 

.  CALL  DEL  3LK<  49 ) 

CALL  BLOCK { 49 ) 

CALL  T  E  X  T  <  F ,  1 ,1000, 10,1  OH DATA  CURVE, TR,F) 


■ 


( I  .  ,0  • 


- 


(IlH.Ih 


'  T  >  \ 


. 


■ 


- 
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CALL  ENDS  LK 
CO  300  I  =  1,NPERD 
SCRENT  =  T  < I ) /XSCALE 

I _ AX(I)  =  ALQG 1C !  SCRENT  ) *  *  ?  0  o _ 

300  AY (  I )  =  V  < I  I  *  1000 
CALL  E RS3K (  1 ) 

CALL  DELBLK(  it  ) 

(CALL  BLOCK ( 10  ) 

C 

C_ DATA  CURVE 


CALL  DLINE(F,F,AX,AY,NPERD,F,F) 


X-AXI S 

CALL  MOVE ( F , 0 ,0  ) _ 

CALL  VECTOR! F  ,800,0, F  ,F  ) 
CALL  ENDBLK 

CALL  DISPLY! 10,1,0, THERE) 
CALL  DISPLY(45, 1,0,0) 

310  CALL  TRANMT 

CALL  E  R  S  B  K ( 52  ) 


A  VALUE  OF  K  MUST  BE  ENTERED  A1  GRID  KEYBOARD  .  .  . 
CALL  C  ALP  HA  (  1  ,  I  X,  I  Y  ,  N IJ  M ,  10, STR ING ,£3207 

_ o  o  .  OTHERWISE  ERROR  MESSAGE  APPEARS _ 

C 

309  CALL  RST3K ! 52  ) 

GO  TO  310 
C 

320  CALL  KR  ( K AY  1  ,R { 1  ),R(2), STRING, NUM , 1 , 2 , £  30  9) 

C _ _ 

C  T 1  IS  ENTERED  AT  TERMINAL  KEYBOARD 

C 

17  WRITE!  o,31) 

81  FORMAT  (  •  ENTES  Tlr) 

READ! 5,93)  T1 

93  FORMAT!  F10. 2  ) _  _ _ 

C 

C 

C _ HI  CALCULATED  AND  PRINTED 

C 

HI  =  SQRT (R(l)*Tl*10#0)/8* 

_ WR I T E ! 6 , 103  )  HI _ _ 

103  FORMAT (6H  Hi  =  ,F10.2) 

HH  *  FP  l=4c  HI  IG 
C _ 

C  FLASHING  *  AT  T  ON  THE  X-AXIS 

C 

_ I  T  1  =  200*AL0G10(  T1  ) /XSCALb _ 

•  CALL  DELBLM3  ) 

CALL  BLOC  K ( 3 ) 

_ _ _  CALL  TEXT<F, IT  1,0,1, 1H*,TR,TR) 


C 

C 

c 

c 

c 
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CALL  ENDBLK 

CALL  D I S P LY ( 3  ,1,0,IHERE) 

C 

c _ 

C  #  If.  if.  if.  if.  #  X'  ❖  £  X  $  $  #  X  Xf  sjs  Xf  Xf  3}c  #  #  XfXr  #  #  ##  #  $  sjcsjt#  sjt  #  #  sjc  sjc  $  $j{s  ^  #  #  ##  #  ijeaj:  sj:#  #  ajc 

0  3}:  jjt  5^  J{!  ijs  ^  ;Jc  3*t  ?Jc  J$:  jJ;  jj:  Xf.  j>;  ^  ;£  jj-  £  >,„•  j|j  ,{;  ^  jj,  jJ-  ;v  }<{  5)4  }*{  #  3^  jjf  ijt  3^  3^3,  3^  3^  jjt  s£  if.  3^  3jf  if  if  3}C  jfcjJl  X-  X  '  Xf 

c _ 

C  WE  NOW  BEGIN  CALCULATION  GF  THE  TWO-LAYER  APPROXIMATION 
C  WHICH  DEPENDS  CN  THE  OBSERVED  VALUE  GF  Tl,  THE  ”R IGHT- 

C  MOST  ZERO-CROSSING  OF  THE  TIME  AXIS  BY  THE  FUNCTION  V. _ 

C  NOW  THAT  WE  HAVE  THE  VALUES  FOR  THE  MODEL,  WE  SHALL 

C  CALCULATE  THE  VALUES  OF  THE  FIRST  APPROXIMATION, 

C  __  WHICH  IS  ESSENTIALLY  A  2-LAYER  CURVE.  THIS  CURVE 
C  SHOULD  FIT  THE  MODEL  CURVE  UP  TO  A  PERIOD  T2. 

C  THE  CALCULATIONS  ARE  STRAIGHTFORWARD  AND  EXACTLY  AS 

C _ LAID  OUT  IN  THE  MGZESQN  THESIS, _ 

C 

C 

C 

C 

1000  DO  1300  J  =  1 , N  PE RD 

_ ARG(J)  =  -HF/SQRT(1Q.*T ( J)  ) _ 

FN(J)  =  EXP< ARG ( J ) )*CCS ( ARGI J ) ) 

FNC(J)  =  EXP( /!RG<  J)  )*C*SIN(ARG(J)  ) 

_  V1CC { J  )  *  KAY1* ( FNC ( J )  +  F 1 4  ( J  )  > 

VT  (  J  )  =  K'AY1*FN(J) 

VI (J)  =  VT( J  ) 

_ V  1C ( J  )  =  VI ( J  ) _ 

R  1  <  J )  =  V  1  ( J  )  -  V(J) 

R1C(J)  =  VIC ( J )  -  V(J) 

1300  CONTINUE 

FRAC  =  ( l.-KAYi) /( 1.+KAY1 ) 


C 

C _ DISPLAY  FIRST  AP  PROX  I  MAT  I  ON 

C 


1800  CALL  DEL BLK ( 49 ) 
j _ _ CALL  BLOCK (49) 

CALL  TEXT ( F , I ,1000,  l  9 , 1  9  H  F  I  R  S  T  A PPR OX  I M AT  I  ON  ,  T R , F ) 

CALL  ENDBLK 

_ CALL  GRID! 1 , T  2 , K A  Y 1 , K  A  Y  2 , VIC, T , R ( 1 )  ,R( 2)  ,R( 3 ) ,uN,NPERD, IHERE, 

2X SC  ALE, 015, Cl  00 0,81900, 820 00 ) 


C 

C  _ T  AND  H  CALCULATED  AND  DISPLAYED 

C 

1900  H2  =  l./FRAC*(SQRT< 10.*R(l)*T2)/8.  -Hi) 

_ WR  I  T  E  (  6,  104)  T  2 ,  H2 _ _ _ 

CALL  TH(T2»H2»IHERF,X$CALE) 

GO  TO  1800 

134  FORMAT (6H  T^  -  . F10«2 , IQXt 6H  H2  =  ,F10.2) 


C 

C 

c 

c 

c 

c 


NEXT,  WE  CALCULATE  THE  THREE-LAYER  APPROXIMATION,  WHICH 
IS  DEPENDENT  'UPON  THE  POINT  AT  WHICH  THE  FIRST  APPROX- 


(  U3H1 

(t  >TV 


i .  u  a  m 


■ 


. 


IV  G  COMPILER 


GRAFIC 


00-17-71 


10:45.39 


PAGL  0006 


C 

c. 

C 

c 


IMATIGN  LEAVES  THE  V-CUkVE,  T2. 


C 

2000 


on  1330  J=1,NPERD 
A R G 1  {  J )  =  A R G ( J ) *FRAC*H2/Hl 
F N 1  {  J )  =  EXP ( ARG 1 ( J ) )  *CGS'lA RG1 ( J  )  ) 
FN1C(J)  =  EXP(ARG1( J ) )*C*SIN(ARG1< J) ) 
VTT(J)  =  KAY2*FN 1 ( J ) 


1330 


V2(J)  =  V1C(J)  +  VTT(J)*FN(J) 

V2CC  (  J  )  =  (V  ICC  {  J  )  +KA  Y2*  (FN1  (  J  )  +FN  1C(J))*(FN(J)  +  FNC(J)  )  )  /  ( 1 .  +  K  AY  1* 
CKAY2*  (FN1  ( -J  )  +FN1C  (  J  )  )  ) 

V2C ( J)  =  RE A  L ( V2CX  ( J )  ) 

R2  { J  )  =  V2 ( J  >  -  V<J) 

R2C { J)  =  V  2C ( J)  -  V  (  J  ) 


CONTINUE 

FRACi  =  (1.0  -  K AY2  ) / (1.0  +  KAY2) 


DISPLAY  SECOND  APPROXIMATION 


CALL  PE  LQ  LK ( 49 ) 


CALL  BLGCK( 49  ) 

CALL  TEXT! F, 1 , 1000, 2O,2OHSEC0ND  APPROX  IM AT  ION , TR , F ) 

CALL  ENDS  LK 

CALL  GRID(2,T3,KAY2,KAY3,V2C,T,R(2),R(3),R(4)  ,0N , NPERD , I  HERE , 
2X SCALE ,£1000, 62000* L2  900, £3000) 


C 

C 


T  AND  H  CALCULATED  AND  DISPLAYED 

^9jG  H3  -  SORT! 10. 0*R( 3) ) /d.* ( SORT ( T 3) -SQkT( T2 ) ) 
WRI T E { 6 ,10 5  )  T3 ,  H3 

105  FORM  AT  (  6Fi  T3  =  ,  F  1 0 . 2 , 1 0  X  ,  6H  H3  =  ,F10.2) 
CALL  TH( T3, H3 , IHERC , X5CALE  ) 


GO  TO  2800 


C 

c 

c 

c 

c _ 

c 

c 

c 

c 

c 

3000 


i  $  $  *  $  £  $  $ $  $  i  $  i  $  .1  $  $  $  $  $ 

if  if  if  >Jc  5;:  if  if  if  if  if  if  if  if  if  if  if  if  if  if  if  if  if  if  if  if  if  if  if  if  if 


$$*$$$$$$$$$$$$$$$$$ 
i-  >|c  if  if  if  if  if  >J;  if  if  if  ;J;  if  if  if  if  ^  if  if  if  if  f  if  yf  if  if 


NOW,  nE  CALCULATE  THE  F JUR-LAYER  APPROXIMATION,  USING  THE 

OBSERVED  VALUE  OF  T3.  WE  WILL  OBTAIN  A  VALUE  FOR  T4, 

THE  POINT  AT  wHICF  THIS  APPROXIMATION  LEAVES  THE  V-CURVE. 


DO  1510  J~l, NPERD _ 

ARG2 ( J  )  =  ARG1( J >*FRAC1*H3/H2 

FN 2 ( J )  =  EXP( ARG2(J  )  )  *CCS ( ARG2 ( J ) ) 

FN2C ( J  )  =  E  XP ( ARG2 ( J )  ) *  S  I N ( APG2 ( J ) ) *C 
VTTT(J)  *  K  A  Y  3  *  F  N  2  (  J ) 

V3 ( J )  =  V2C ( J )  +  VTTT(J)*FN( J)*FN1( J) 

V  3  C  C ( J )  =  (V1CC(J)  +  KAY2*(FNi(J)  +  FN  1C ( J )  )  * ( FN ( J )  +  F  N  C  <  J  )  )  + 
.CKAY3*(  FNi  (  J  )  +  FN1C<  J  )  )  *(  FN<  J  )  +  PNC ( J ) ) *( FN2 ( J )  +  FN2C(J)) 
c+  V  1CC(  J  )*KAY2*KAY3*  (FN2  (  J  )  +  FN  2C  (  J  )  )  )  /  (  1  •  0  ♦  KAY1*KAY2* 

C(FNi(J)  +  FN1C(J))  +  KAY1 *KAY'i* ( FN 1 ( J )  +  FN1C  (  J )  ) * ( FN2 ( J  )  +  FN2C ( J )  ) 


I  ) 


’ 


w  (  l  >  1  ) 


h,  t  riSw 


'  ■.  JbM  *A>  3j?  i  ‘  Hi 


i  i  \  >  1  > 


. 

. 

1  )  (  Cl  > 


*  i  t,  )  ,  )  ; 
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C  +  KAY2* *KAY3* (FN2(  J )  +  FN2CIJ))) 

V3C ( J )  =  REAL ( V3CC ( J ) ) 

R3(J)  =  V  3  <  J  )  -  V ( J ) 

_ R3C  (  J  )  =  V3C  (  J  )  -  V  (  J  ) _ _ _ 

1510  CONTINUE 

FRAC2  =  ( 1 .0-KAY3 ) /{ 1 .0+KAY3 ) 

€ _ _ _  _ 

C  DISPLAY  THIRD  APPRGI XMATI CN 

C 

3800  CALL  DEL BLK  (  45  )  _ _ _ 

CALL  BLOCK (49) 

CALL  TEXT ( F ,  1 , 1000,  1 9,19HTHIRD  APPROX IMAT I C N , T R , F ) 

CALL  EN  3BJL  K 

CALL  GR 10 ( 3,T4, KAY 3 ,KAY4, V3C, T,R ( 3)  ,R( 4) ,R( 5 ) fGN, NPERO, I  HERE , 
2XSCALE, £2000,  £3000, £3900, £4000) 

C _ _ _ 

C  T  AND  H  CALCULATED  AND  DISPLAYED 

C 


3900  R (4 )  -  R (3)  /FRAC2 

H4  =  SORT (R (4 )/R ( 1 )  MU  SORT {  iQ.*R( 1) *T4> /8.  -FRAC*H2  -HI 
C-FRAC*FRAC1^H3 ) 

_ WRITE (6. 110)  T4,H4 _ 

110  FORMATION  T4  =  , F 1 0 . 2 ,  10 X , 6H  H4  =  ,F10.2) 

CALL  TH ( T4,H4,IHERE, XSCALF) 

GO  TO  3300 


C 

c 

c. 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


>>»> . >»» . >>>>>..  >>>». ....  »»> . >»>> 


FINALLY,  WE  SHALL  CALCULATE  THE  FIVE-LAYER  APPROXIMATION, 
USING  THE  OBSERVED  VALUE  CF  T4,  THE  POINT  AT  WHICH 
THE  THIRD  APPROXIMATION  LEFT  THE  V-CURVE. 


4000  DO  1600  J  =  i , N  FE  RD 

_ ARG3 ( J )  =  ARG2( J  )^FRAC2*H4/H3 _ 

F N3 ( J )  =  EXP( ARG3( J  )  ) *CGS ( ARG3( J ) ) 

FN  3C ( J )  =  EXP( ARG3( J ) ) *SIN ( ARG3 ( J ) ) *C 
VTTTT(J)  =  KAY4*FN3(J> 

V4(J)  =  V  3  C ( J )  +  V TTTT < J ) *FN( J)* FN 1 ( J )*FN2< J ) 

V4CC ( J )  =  (V1CC(J)  +  KAY2/KAY1*V1CC ( J ) *( FN1( J)  +  FNIC(J))  + 

_ CKAY3/KAY  1*V1CC<  J  )*(  FN1<  J  )  +  FN  1 C  (  J  )  )  *  (  FN  2  (  J  )  .  +  FN2C(J)L_+_ - 

C  KAY4/KAY 1* V  ICC ( J )* ( F M ( J  )  +  FN1C(J))JMFN2(J)  +  FN2C(J))*(FN2(J) 
CFN3CIJ))  +  KAY2*K  AY3*  V  iCC  (  J  )  *  (  FN2  (  J  )  +  FN2CIJ))  +  V1CC(J)* 

C  ( F N 2  (  J  )  +  FN 2 C  (  J  )  )*(  FN3(  J  )  +  FN3  C  (  J  )  )  *KA  Y2*K  AY4  +  V1CC.U1  ♦KAY3 * 

CKAY4* (FN 3 ( J )  +  FN3C(J))  +  KAY2*KAY3*KAY4/KAY1*V1CC ( J I* (FN1C( J I 
C+FN1 ( J ) )* ( FN3( J )  +  FN3C ( J ) )  ) / ( !•  +  KAY2* KAY  1* ( FN1 ( J )  +  FNiC(J)) 

C  +  KAY 1  AY3*  (  FN  1  (  J  )  +  FN1 C  (  J  )  )  *  (  FN2  (  0  )  +  FN2C(J))  +  KAY2*KAY3* _ 

•  C ( FN  2 ( J  )  +  F N2C ( J  )  )  +  K A Y 1*K AY4* ( FN 1 ( J )  +  FN1C ( J ) ) * ( FN2 ( J )  + 

CFN2 C ( J )  ) * ( FN3 ( J  )  +  FN3C(J))  +  KAY2*KAY4* ( FN2 ( J )  +  FN2C(J))* 
C(FN3( J)  +  FN3C(J)  )  +  K AY3*KA Y4* ( F N3 ( J )  +  FN3C(J)  )  +  K A Y 1*KA Y 2# 


♦'.fir  r  +; 


it  t  i *  »  )  r > 


,  (  }  »  <  *  .  <  *  r  «  *  •  1  »  .  >  I  -l.»' 


■ 


r 


' 


->■ 


t'll'a  > 


' 

i  .  >  \  1  • )  1 


. 


XY  1  ♦  lit))  *  ♦  IDS 

)  )  (LI  +  f  L  )  I  *  (  v  > 

1  <  (  *  )  '  1  D  1  ♦  l,  • 
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CKAY3*KAY4*  (  FN  1  (  J)  +  F N1 C ( J  )  ) * ( FN3 ( J  )  +  FN3C(J))) 
V4C ( J  )  =  REAL (V4CC{ J )  > 

R4{ j)  =  V 4 ( J  )  -  V( J  ) 

R4C ( J  )  =  V4CCJ)  -  V ( J  ) _ 

1600  CONTINUE 


C 

C 

C 


DISPLAY  FOURTH  APPROXIMATION 


CALL  DELBLM45) 

_ CALL  BLOCK (49) _ 

CALL  TEXT { F , 1 , 1000, 20 ,20HF0URTH  APPROX IMAT  ION , TR , F ) 

CALL  ENDBLK 

CALL  GRID (4 ,T4, KAY4, KAY 5, V4C, T, R (4) *R{ 5 ) ,R(6 ) , ON ,NPERU, I  HERE , 
2X SCALE ,£2000,63000, 6995, £999) 

C  JOB  TERMINATES  NORMALLY  IF  ANY  ATTEMPT  TO  CONTINUE 

999  RETURN _ _  __ 

END 


EMORY  REQUIREMENTS  008274  BYTES 
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SUBROUTINE  SCREEN 
C 

C  SUBROUTINE  SCREEN  INITIALIZES  THE  GRID  DISPLAY 

C 

LOG  I C  AL* 1  TR/. TRUE. /, F/. FALSE./ 

C 

C  BLOCK  I  DISPLAYS  THE  APPROXIMATION  CURVE 

C 

CALL  BLOCK ( 1 ) 

_ CALL  T6XT!F, 1 ,1,1,1HA,F,F) _ 

CALL  ENDBLK 
C 

C  BLOCK  2  DISPLAYS  VALUES  OF  KAY,R,T,H 

C 

CALL  BLOCK ( 2 ) 

_ CALL  TEXT(F,2  ,1,1»1HA,F,F) _ 

CALL  ENDBLK 
C 

C  BLOCK  3  DISPLAYS  A  FLASHING  *  ON  THE  X-AXIS  AT  T 

C 

CALL  BLOCK ( 3 ) 

_ CALL  T£XT!F,1 ,100,1, 1H*,F,F) _ 

CALL  ENCBLK 
C 

C  BLOCK  10  DISPLAYS  DATA  CURVE  AND  X-AXIS 

C 

CALL  BLOCK! 10  ) 

_ CALL  TEXTIF, 1  ,1, 1 , 1HA ,F,F) _ 

CALL  ENCBLK 
C 

C  BLOCK  49  INDICATES  I  HE  APPROX IMATH  i  DISPLAYED 

C 

CALL  BLOCK ( 49 ) 

_ CALL  TEXT ( F,  1  , 1000*20  ,2QHZERQTH  APPROX  I M AT  I  ON , T R , F  ) _ 

CALL  ENCBLK 
C 

C  BLOCK  50  DISPLAYS  INDICATERS  AND  INSTRUCTIONS 

C 

CALL  BLOCK ( 50 ) 

_ _ CALL  TEXT(F,-809,-40,47,47H  KAY _ R _ I _ 

2  H,TR,F) 

CALL  TEXT (F ,0  ,0 ,5 ,5HCLEAR, TR,F ) 

CALL  T  EXT ! F  ,  C  ,-40 , 4  , .4 H BACK  ,  TR  ,  F  ) 

CALL  T  EXT ( Ft  0  80  ♦ 9 , 9HHAR0  COPY,TR,F) 

CALL  TEXT ( F  , 0  ,-l 20 , 8 , 8HCGNT I NUE , TR , F ) 

_ _ CALL  TEXTIF, 0 ,-160, 12, 12HBLINK  APPROX,  TR,F) _ _ 

CALL  ENDBLK 
C 

C  BLOCK  51  ALLC'aS  USER  TG  RESET  BLOCK  5U  IF  IT  IS  CLEARED 

C 

CALL  BLOCK!  51  ) 

_ CALL  TEXT (FfO  5,5HR E SET fTR>F) _ 

' CALL  ENCBLK 
C 

C  BLOCK  52  INDICATES  ERROR  BY  USER 


.  va 


LI*  L  JJAi 


(!)  >.i  IJ 


•  • 


,  r  U  r  '  1  - 


.  1  .  -t  e  I T  T  1.3 


tt  JLT.1 JU5IJLL*  f  JJA3  „ 


. 

■a 


- 

. 
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C 

CALL  BLOCK  < 52  ) 

CALL  TEXT ( F , 397 , 5 12, 17 , 1 7HERR0R  -  TRY  AGAIN, TR,TR) 

I _ CALL  ENJBLK _ _  _ 

C 

CALL  D  I  S  PL Y ( 1  ,1 , 0,0) 

CALL  C  i  i>P  l  Y  {  2 ,1,0,0) 

CALL  DISPLYI 3, 1,0,0 ) 

CALL  DISPLYI  1 C  ,  1 , C , 0 ) 

I _ CALL  OISPLY  (50ti,810«  1000) _ 

CALL  DISPLYt 51,1 ,810,1000) 

CALL  DISPLY <52, 1,0,0) 

CALL  ERSBK  <  51 ) 

CALL  E  RS  BK ( 52 ) 

RETURN 

END 


MEMORY  REQUIREMENTS  00055A  BYTES 


■ 

' 
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SUBROUTINE  GR ID( I  I , T I , KA Y , KA YN , V , T , RP , R , RN ,QN , NPERD , I  HERE , XSC AL E , 

2  *,*,*,* } 


C 

C 

c 

c 


c 

c 

c 


c 

c 


c 

c 

c 

c 


c 

c 


c 

c 

c 

c 

c 


c 

c 

c 

c 


SUBROUTINE  GRID  DECODES  AND  PERFORMS _ 

INSTRUCTIONS  FROM  THF  GRID  USER 

INTEGER*?  AX< 100) , AY l IOC) , STRING! 10) 

INTEGER  BLK(  8  )  , STATUS, KEY 
REAL  KAY,TI ,T ( 100 ) ,V< 10Q) 

LOG  I  CAL*  l — QN  ,  ERR  ,  F  /.  FAL  SE  ./,  Tk  /.  TR  U  L.  /«  F  I  A  SH  /.)-  A  t_  SC  «  / 
APPROXIMATION  CURVE  IS  DISPLAYED 


WRITE (6,103)  II 

103  FORMAT ( 15H  APPROX IMAT ION  ,11, 14H  NOW  DISPLAYED) 

II  I  =  I  I  +  1 _ _ 

CALL  DISPLY (49,1,0,0) 

DO  5  I  =  1,  NPERD 

SC RENT  =  T (  I  )  /XSCALE 

AX (  I  )  =  AL0G1 0( SCRENT)*20G 

5  AY (  I  )  = V (  I  )*10C0 

6  CALL  DEL BLK  (  i  ) _ _ 

CALL  BLOCK!  1  ) 

CALL  DLINECF ,F, AX, AY, NPERD, FLASH, F) 

CALL  E NQ p-  LK 

CALL  0 1 SPLY (1  *  1 , 0 , 1 HERE) 

CONTROL  TRANSFERRED  TG  GRID  USER _ „ 

10  CALL  TRANMT 

_  ERROR  INDICATION  ERASED  IF  NECESSARY 

CALL  ERSBM  52  ) 

GO  TO  200  IF  LIGhT  PEN  USED _ 

CALL  DLPEN( 1 , IX, I  V, I  TYPE, ID, BLK, £200) 

GO  TO  500  IF  KEYBOARD  .  .  •  OTHERWISE  ERROR 

Call  DALPHAC 1vIX« IY,NUM, 10,STRI <G,£50  ) 

9  CALL  RSTBM52) 

GO  TO  U 

LIGHT  PEN 


GO  TG  300  IF  USER  HAS  PICKED  A  POINT  ON  THE  GRAPH 
.  •  .  THAT  IS,  WISHES  TO  CHANGE  T 

200  IF( IX.LT.810)  GO  TO  SCO 

CLEAR  OR  RESET  INSTRUCTIONS 
I F ( I Y.GT.93Q  )  GO  TG  21a 

GO  BACK  TO  PREVIOUS  APPROXIMATION  OR  DATA  CURVE 
I  F  (  IY.GT.9V0)  RE  TURN  l _ _ _ 


C 

C 


HARD  COPY 

IF ( IY  .GT .900 )  GO  TO  2  40 


UdOU  TCtl 


i  (  , 


, 

- 
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C 

C  CONTINUE  ON  TC  NEXT  APPROXIMATION 

I F ( IY.GT.860)  GC  TC  250 

C _ _ _ _ _ _ 

C  FLASH  APPROXIMATION  CURVE 

IF ( I Y.GT. 820 )  GO  TO  2b0 


C  O.EAR/RESET 

C 

210  I  F(  ON )  GC  TC  215 


CALL  ERSBM51) 
CALL  RST BK ( 50 ) 
CALL  RSTBM 2 ) 
ON  =  .TRUE. 

GO  TO  10 

215  CALL  ER SBK ( 5 C  ) 

CALL  RSTBM  51  ) 
CALL  ERSBK ( 2 ) 
ON  =  .FALSC. 

GG  TO  10  ' 

C 

C  HARD  COPY 


240  CALL  ERSBK ( 50  ) 

CALL  ERSBK( 2 ) 

CALL  SNAP (0.0 15) 

CALL  RST8K( 50 ) 

_ CALL  RST  3K  (  2  ) _ _ 

GO  TO  10 

C 

C  CONTINUE 

C 

250  CALL  DALPHA(2, IX, IY,NUM, 10, STRING, £255) 

_ GO  TO  9 _ _ _ 

C  KAY  FOR  NEXT  APPROXIMATION  MUST  HAVE  BEEN  ENTERED 

C  AT  GRID  KEYBOARD 

25  5 _ CALL  KRIKAYNt R«RN« STRING »NUM, 11+1,111  +  1,09) 

RETURN4 

C 

C _ MAKE  THE  APPROXIMATION  CURVE  BLINK  OR  STOP  BLINKING 

C 

2bU  FLASH  =  .NOT. FLASH 
J _  GO  T  Q_6 _ 

C 

C  NEW  VALUE  FOR  T 

C _ _  _ 

C  •  CHECK  IF  USER  H 4 S  TYPED  IN  A  VALUE  FOR  T  .  .  . 

300  CALL  0 ALP HA (2, IX, IY,NUM, 10, STRING, 6350) 

C  ...  IF  NOT,  THEN  USE  POSITION  UF  LIGHT  PEN 


;  .  7  i  -%.£  it»j»  *8  bf  i  ■'  <  *  .  At  ■ '  '«:.to 
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360  T I  =  XSCALE*10**{  IX/200.  ) 

RETURN  3 

C  ...  IF  SC,  THEN  CONVERT  AND  USE  VALUE  FROM  GRID  KEYBOARD 

1~.  350  I  F(Q.  EQ.KODC  NU.V.2  )  >  NQM  =  2*NUM _ 

IF  ( 1 • EQ. MOD ( NUM, 2 ) )  NUM  =  2*NUM-1 
T I  =  CHRFLT (STRING, NUM, ERR ) 

IF (ERR)  GO  TO  9 
RETURNS 
C 

C _ KAY  ENTERED  AT  KEYBOARD _ 

C 

500  CALL  KR( KAY, RP, R, STRING, NUM , I  I , I  I  I , 09 > 

_ RE T UR  N  2 

END 

gMOEV  RFQUI  RFM  FNTS  OOP  AC2  BYTES 
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SL8R0UT INE  KR ( KA Y,RP, R, STRING »NUM»II » I II ,*) 

C 

C  KR  DECODES  K AY t C A LCUL ATES  R  AND  DISPLAYS  THEIR  VALUES 

C _ 

LOGICAL* 1  TR/.  TRUE, /, F/. FALSE. / »  ERR 
I NTEGER*2  STRING! 10  ) 

REAL  KAY  » R 

NUM  =  2 * NUM  -  MOD ( NUM t  2 > 

500  KAY  =  CHRFLT (STRING, NUM, ERR) 

_ _ IF(ERR)  RETURN! _ _ 

FRAC  =  (1.0  -  KAY )  / (1.0  +  KAY) 

R  =  RP/FRAC**2 

WRITE (6, 102)  I  I,  KAY, I  II , R 

102  FORMAT  1 4  H  KAY,I1,3H  =  ,  F10.3 , 3.0X  ,  3H  R  (  ,  1 1  ,4H  )  =  ,ii  (,  ) 

CALL  DELBLK ( 2 ) 

I _ CALL  8LQCK(  2  ) _ _ _ _ 

CALL  TEXT(F,1  ,930,6,STRING,TR,F ) 

CALL  FLTCHR ( R, STRING, 8,2) 

CALL  TEXT(F,192,930t 8 T STRING, TR,F) 

CALL  ENOBLK 

CALL  DISPLY(2, 1,0,0) 

K _ RETURN _ ' _ . _ 

END 
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SUBROUTINE  TMT  ,H,  I  HERE  ,  X  SCAL  E  ) 

TH  DISPLAYS  T  AND  H  AND  FLASHING  *  AT  T  ON  THE  X-AXIS 


LOG  I  CAL*  1  TR/.TRUE./, F/. FALSE./ 
INTEGERS  STRING!  10  ) 

CALL  DEL BLK ( 2 ) 

CALL  BLOCK! 2) 

CALL  FLTCHR!T,STRING,9,2 ) 

CALL  TEXT(F.41b.930t9,STRING,TR.F> 

CALL  FLTCFR(H,STRING,8,2) 

CALL  TEXT(F ,672 ,930,8 , ST RING, TR,F ) 

CALL  ENDBLK 

CALL  DISPLY(2, 1,0,0) 

CALL  DE  LB  LK ( 3 ) 

CALL  BLOCK ( 3  ) _ 

ITX  =  200  *  ALGGICIT/XSCALE) 

CALL  TEXT!F,ITX,0,1,1H*,TR,TR> 

CALL  ENDBLK 

CALL  DISPLY! 3,1,0 , IHERE ) 

RETURN 

END 
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